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1.0 SUMMARY 


This report presents a panel-type influence coefficient method which solves the linear integral equations 
of unsteady, small disturbance, subsonic, potential flows. The method is applicable to flows about 
arbitrary, nonplanar bodies having flow boundary surfaces undergoing small amplitude harmonic 
oscillations. The flow boundary surfaces are paneled (i.e., subdivided into small quadrilateral 
surface segments) and polynomial distribution splines of unsteady sources and doublets are 
defined for each panel. The method is an outgrowth of a recently developed steady flow panel 
method and employs the linear source and quadratic doublet splines of the steady flow method. 

Since the flow boundary surfaces are undergoing small amplitude harmonic motion, they have 
mean steady locations relative to which their unsteady motion is described. As a consequence, 
the unsteady flow problem is formulated as a small disturbance to a mean steady flow, and the 
mean steady flow is the flow about the arbitrary, nonplanar boundary surfaces in their steady mean 
locations. This formulation of the unsteady flow problem follows from a Taylor series expansion of 
the unsteady flow boundary conditions about the steady mean flow boundary conditions. The 
result is an unsteady flow problem which depends on an underlying steady mean component of flow 
and which contains boundary conditions evaluated at the steady mean locations of the boundary 
surfaces. The boundary conditions are such that the unsteady flow analyst can incorporate important 
features into the flow model which are not available in the usual theory of linearized flow. 

One of the important features of the formulation is the wake model; the wake can be an 
arbitrarily shaped surface along which the unsteady vorticity can be convected arbitrarily. This 
feature is useful when modeling the unsteady flow about complex configurations, for example, 
those consisting of wing-body-tail combinations. When the usual wake model, consisting of free 
stream convection of the wake vorticity, is employed for such a configuration, the tail cannot 
be correctly positioned relative to vorticity in the wing wake. Using the present panel method, 
if the analyst can determine the appropriate steady mean wing wake surface streamlines which 
follow the contour of the steady mean body surface, then he can solve the problem with the 
unsteady vorticity in the wing wake convected along those streamlines using the steady mean 
convection speed. Using this wing wake model rather than that of the usual linear theory, the 
unsteady wing wake vorticity location and phase relative to the tail more closely represents 
the physics of the problem. 

A second important feature of the formulation stems from the geometric location of the unsteady 
source and doublet distributions; they can be located on the actual surfaces of thick bodies in 
their steady mean locations. The usual linear theory places these distributions on mean defining 
surfaces which are restricted to being cylindrical with the generator of these cylindrical surfaces 
aligned with the undisturbed freestream. This restriction is a consequence of having expanded 
the boundary conditions about the freestream. In the present panel method, the boundary con- 
ditions are expanded about the steady mean component of flow; hence, the sources and doublets 
are distributed on the boundary surfaces in their steady mean locations. Since the unsteady flow 
disturbances are emitted by the unsteady source and doublet distributions, the panel method 
allows the unsteady disturbances to be emitted from locations which are significantly closer to 
their actual physical locations. Even so, the panel method still suffers from another defect of the 
usual linear theory, namely, that disturbances propagate with the freestream speed of sound. 

Regardless of this error, however, the panel method can circumvent part of the error in the com- 
putation of the cause-and-effect phase relation of unsteady disturbances emitted from the 
boundary surfaces of thick bodies. 



Only the second of the aforementioned features of the panel method is demonstrated in the report, 
namely, that thick aerodynamic surfaces can be paneled on their actual, mean steady locations. The 
first feature, which allows a more accurate wake model, is not demonstrated because the steady flow 
computational capability required for computing the mean steady flow wake location and vorticity 
convection was not available. The unsteady flow about a wing-body-tail-nacelle configuration has 
been computed and the wake surface was deformed so that it follows the body contours; however, the 
wake deformation is the result of an engineering judgement and was not a result of an actual steady 
flow computation. 

The primary objective of the work was to demonstrate the capability for modeling unsteady flow 
about realistic aircraft geometries using paneling on the actual surfaces of thick bodies. That objective 
was achieved. 
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2.0 INTRODUCTION 


This report describes a panel method for evaluating the unsteady, subsonic flow about a solid body 
which is of arbitrary shape and which is undergoing harmonic oscillation in a uniform freestream. 

The method is an outgrowth of the steady flow panel method described by Reference 1 , “A Higher 
Order Panel Method for Linearized Supersonic Flow.” The method has been implemented as a 
computer program by modifying the computer program of ref. 1 , but only the subsonic capability 
of that computer program has been extended to unsteady flow. The modification is such that, if 
the unsteady flow program is executed for a harmonic frequency of zero, the two computer programs 
yield identical results. The validity of the unsteady flow panel method has been tested by comparing 
computed results with wind tunnel test results and with results computed by other subsonic, unsteady 
flow prediction methods. Those comparisons are shown in sec. 4 of this report. Sec. 5 contains 
a supplement to the user’s manual for the method of ref, 1 , i.e., “User’s Manual, Subsonic/Supersonic 
Advanced Panel Pilot Code,” ref. 2. From the point of view of the program user, the two programs are 
so similar that a separate user’s manual for the unsteady panel method is neither necessary or desirable. 

The unsteady flow panel method is described in sec. 3 with references to appendices containing the 
details of its derivation. The description is divided into five parts, viz., sec. 3.1 - sec. 3.5. Sec. 3.1 
introduces coordinate systems and the method used to describe unsteady surface motion. Sec. 3.2 
describes the unsteady flow boundary value problem. That boundary value problem is derived in 
App. B as a linear unsteady flow theory approximating the general theory of unsteady flow 
appearing in App. A.f Sec. 3.3 describes the integral equation, which represents a general 
solution to the unsteady flow boundary value problem. The integral equation is derived 
from Helmholtz’s theorem in app. C. 

Sec. 3.4, sec. 3.5, and app. D, together, describe the panel method of reducing the integral 
equation to a determinate system of linear, complex algebraic equations. The solution to this 
system of equations provides a solution to the unsteady flow problem. 

In essence the unsteady panel method is a method for constructing approximate solutions to Helmholtz’s 
equation, viz.. 


V2 0* + (iJ2)2 0* = 0 

where 0* is the complex amplitude of a potential and F2 is a constant referred to as the Helmholtz 
coefficient. 


The solution is constructed in a domain, V , enclosed by a surface, S , by means of Helmholtz’s 
theorem, ref. 3, viz.. 


J_ 

47T 



a /e' 


•i^2R^ 


ds - 0* (P) 


where a *(0) and ju*(Q) represent, respectively, the complex amplitudes of sources and doublets 
distributed on the surface, 2 , and 


r=Ir'-rI 

t This development is significant in that it leads to the so-called “nonlinear boundary 
conditions.” 
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where R is the position of a field point P while R' is the position of a point Q on the surface of 
integration 2 . Also, 

)=ft*v( ) 

on 

where n is a unit vector which is normal to S and which is directed from a point Q on 2 and into 
the volume enclosed by 2 , see fig. 1 . 



Figure 1. — Surface of Integration and Domain of Dependence 


The theory of unsteady flow underlying the panel method is based on the linear, small disturbance 
flow equation. This flow equation, as shown in sec. 3.2, can be related to Helmholtz’s equation by 
changes in the dependent and independent variables. Helmholtz’s theorem, therefore, provides a 
general solution to the unsteady flow equation. A particular solution is obtained by requiring 
to satisfy boundary conditions which are usually specified on the surface of integration 2 . 

As noted above, the panel method reduces the integral contained in Helmholtz’s theorem to a 
linear, complex algebraic expression. As shown in sec. 3.3, for unsteady flow about an aircraft in 
an unbounded atmosphere, Helmholtz’s theorem reduces to a surface integral over the aerodynamic 
and wake surfaces of the aircraft. The quantities 

1 0~i^2R 19 / e"i^R\ 

— and — — I ~ ) 

Air R 47t 9n\R / 
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appearing in Helmholtz’s theorem are kernel functions describing the complex amplitudes of the 
potentials induced at P due to, respectively, a unit source and a unit doublet located at Q . As 
noted above, the complex quantities 

a* (Q) and (Q) 

multiplying those kernel functions are the distribution strengths of the sources and doublets 
distributed on 2 . These distribution strengths are the unknowns of the flow problem and they 
are determined by requiring 0*(P) to satisfy the boundary conditions of the flow problem. The panel 
method is applied by subdividing the surface of integration into quadrilateral panels, fig. 2. The 
source and doublet distributions are approximated on each panel by a truncated power series, e.g., 

O* (Q) = O* + ^ + T] 

where ^ , t? are the coordinates of a surface point Q and where a^* , a^* , are the complex 
coefficients of the terms of the truncated power series. The surface integral is now evaluated for 
each panel to obtain an expression for 0*(P). This expression is a complex algebraic equation 
in terms of the power series coefficients. The power series coefficients are related by a least 
squares fit to the values of a* and p* at discrete points on the surface. These values are 
collectively called the singularity parameters and they are the unknowns of the problem. 

The flow problem can be expressed symbolically in terms of the singularity parameters as 


<P* (P) = 2 Aj* (P) Xj* 

j 



Figure 2. — Typical Configuration Surface Pane! Arrangement 
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where Xj* denotes the jth singularity parameter. When the values of this expression or of the 
gradient of this expression are specified by the boundary conditions at a sufficient number of points, 
a determinate system may be obtained, for example, 

j 

and that system may be solved for the unknown singularity parameters. The flow characteristics 
which constitute a solution to the unsteady flow problem are then computed making use of the 
computed values of the singularity parameters. 

The panel approximations employed in the present, unsteady flow panel method are identical with 
those of the steady flow panel method for ref. 1 . In general, source distributions vary linearly and 
doublet distributions vary quadratically on each panel; but, depending on the distance of the field 
point from the panel, the linear and quadratic distributions are defined over: (1) eight triangular 
subpanels of the quadrilateral panel, (2) two triangular subpanels of the quadrilateral panel, or (3) the 
quadrilateral panel itself. In any case the aerodynamic influence coefficients, , are evaluated by 
carrying out the integrations, appearing in Helmholtz’s theorem. These integrations are carried 
out independently over each panel region for which the source and doublet distribution ap- 
proximations are defined. This operation leads to values of potential and potential gradiant at field 
points where boundary conditions are applied. The resulting integrals relate these field point values 
to unit values of the coefficients in the power series approximations to the source and doublet 
distributions; they, in turn, are related by the least squares fit to the singularity parameters. 

The integrals involved in evaluating the influence coefficients are decomposed in such a manner, 
app. D, that every integral which must be evaluated is related to two fundamental integrals. The 
integrands of these fundamental integrals are approximated by a series of constrained Tschebycheff 
polynomials. The result is a sequence of integrals which are integrable in closed form. The number of 
series terms required to acliieve a specified accuracy for this approximation depends on the magnitude 
of the Helmholtz coefficient, , and on the panel size; an error analysis of the approximation 
appears in app. D. The series approximations are such that the integrands and the limits of 
integration appearing in the series terms are independent of the Helmholtz coefficient; thus, the 
formulation, like that of ref. 4, would allow the integrals to be evaluated, saved, and recombined 
sequentially to produce aerodynamic influence coefficients for a sequence of Helmholtz coefficient 
values (i.e., a sequence of harmonic frequency values). The computer code, however, was not 
structured to allow this type of sequential computation. 



6 


3.0 UNSTEADY FLOW PANEL METHOD 


3.1 COORDINATE SYSTEMS AND SURFACE MOTION 

The unsteady tlow boundary value problem and the panel method for its solution are described in 
terms of two coordinate systems: 

(1) Compressibility coordinate system - x, y, z 

(2) Local panel coordinate system - ^ ^ 

The compressibility coordinate system is a rectangular, Cartesian system which is an inertial reference 
trame translating with a steady velocity in the negative x direction relative to the undisturbed fluid, 
fig. 3; thus, there appears to be a uniform freestream in the positive x direction. 

As noted in sec. 2 and illustrated by fig. 2, the aerodynamic surfaces and the wake surfaces are 
subdivided into a large number of segments called panels. Each segment of surface represented by a panel 
is approximated by one or more flat panels. A local panel axis system is defined for each flat panel. 

This coordinate system is related to the compressibility coordinates by a transformation such 
that the r? coordinate plane is paraded with the plane of the flat panel and such that the f axis 
is normal in relation to it, fig. 4. The coordinate transformation, relating the local panel' and 
compressibility systems, is described in app.C. 


z 




The above notation (i.e., x , y , z ) for denoting the coordinate lines of the compressibility system 
is different from that of ref. 2, “User’s Manual, Subsonic/Supersonic Advanced Panel Pilot Code.’’ 

Since ref. 2 is used in conjunction with sec. 5 of this report as a user’s manual for the unsteady 
panel method code, this difference in notation should be carefully recognized. Ref. 2 employs x , y , z 
to denote the coordinate lines of the global coordinate system which is used as a basis for describing 
body surfaces and flow properties. The coordinate lines of the compressibility system are denoted 
by Xc V yc , Zc . In this report the global coordinate system is mentioned only in sec. 5. 

The analytical statement of the flow boundary conditions requires an analytical description 
of the aerodynamic and wake surfaces. This description is expressed in terms of surface points defined by 
surface coordinates which are the local panel coordinates introduced above. The surface coordinates 
appear as the parameters in the Gaussian description of the surface (ref. 5, p. 183), viz., 


R = R (g,7?,t). 


If Tj and t are held fixed while ^ is varied, then the position vector, R , traces out a spatial 
curve on the surface as the surface appears at the instant of time t . This curve is a surface coordinate 
line corresponding with the | parameter. Surface coordinate lines corresponding with the tj parameter 
are constructed by holding ^ and t fixed, and varying r\ . A surface coordinate system is seen to 
be formed by ^ and t? and a surface point is defined by a pair of values for | , tj . 

Unsteady motions of the surfaces are described by an unsteady displacement of the surfaces from 
steady mean locations, fig. 5. Let Q denote an arbitrary surface point of the steady mean location 



Figure 5. — Surface Displacement 
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of a surface. This point has the steady position R ' relative to the compressibility coordinate system. 
Unsteady surface mo^n is introduced by allowing this surface point to have the unsteady (i.e., time 
dependent) position R . The surface point therefore undergoes the following time dependent 
displacement: 

D = R-R'. (1) 


The velocity of the surface point, viz., 



( 2 ) 


where g , rj are held fixed, is then used to compute the speed of displacement of the surface in the 
direction of the surface normal, i.e.. 


A 

Un = n 


ao 

at 


(3) 


where n is the unit vector normal to the moving surface at the surface point. 


The unsteady flow theory used in the panel method is an approximate theory which, as shown in 
app. B, follows by assuming that the displacement, D , is small everywhere on the surfaces and by 
assuming that the harmonic motion, i.e.. 


- — = icoD*ei^t 
Uo 9t 


where co = ’ 


is such that the frequency of the motion satisfies 

|dJD*|<G(|D* I); (4) 

hence. 


From the assumed smallness 

0 ^-(vX d) 

n ris + 0 X ns 

A 9D 
u„«^ n • — 


cD<0(l). 

of D the following approximations are taken to be valid (cf. app. B): 
is a first order approximation to surface rotation, 

is the unit vector normal to the moving surface, and 
is the speed of displacement of the moving surface. 
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3.2 UNSTEADY FLOW BOUNDARY VALUE PROBLEM 


The flow is assumed to be compressible, inviscid, irrotational, and isentropic. Under these assumptions 
the motion of the fluid may be described in terms of a velocity potential satisfying the weifknown, 
nonlinear flow equation shown as equation (A. 17) of app. A. In the approximate theory, however, 
disturbances to the uniform freestream are assumed to be sufficiently small that the nonlinear flow 
equation may be approximated by the small disturbance flow equation, viz.. 


V 




Dt^ 


( 5 ) 


where 


_D ^ 

Dt ~ at ® 3x 


and <j) is the disturbance velocity potential which is related to the flow velocity as follows; 

V^ = Uo+v<A. 


( 6 ) 


The panel method is predicated on the separation of the velocity potential into two components, viz,, 


where the velocity potential <j)^ is independent of time and describes the flow if the speed of 
displacement, Un , of every surface point vanishes for a sufficiently long period of time. The 
velocity potential is that of the unsteady component of flow arising as a consequence of 
nonzero Uj^ . 

The steady velocity potential is required to satisfy the steady form of eq. (5), viz., 

(' ■ Wyy ^ 

and the boundary conditions of the approximate theory of steady flow described by ref. 1 . Those 
boundary conditions are as follows: 


Wg • Ug = O onSgand Wg (9) 

where (ref. 1 , sec 3.10) Sg and Wg denote, respectively, the aerodynamic and wake surfaces. 

The aerodynamic and wake surfaces, injthe present application, coincide with the steady mean 
locations of those surfaces. The vector represents the mass flux of the steady component of 
flow which is given the following first order approximation (ref. 1 , sec. 3.1): 


where, in turn, 


pQ^s= Wx) 

and the quantity is the perturbation mass flux of the steady component of flow. 

In addition to this boundary condition, across wake surfaces the pressure coefficient is required to be 
continuous, i.e., 

[Cp]=0onW^. (12) 


The boundary condition shown as eq. (12) is applied using the following second order approximation 
for the pressure coefficient: 







(13) 


This approximation is chosen because (as shown by sec. 3.3, ref.l) it satisfies Euler’s equations 
so that the momentum integral (ref. 1 , eq. 24) is valid; as a consequence, the aerodynamic force 
is computed as follows: 




(14a) 


when the pressure coefficient is approximated by eq. (13). Similarly, the moment of the aerodynamic 
pressure is computed as 


Ms=-P 




CpR^Xfi^ds 


(14b) 


where Rg describes position on the mean steady surface and n^ is the unit normal. 

The unsteady potential, must satisfy eq. (5) and the following boundary conditions 
derived in app. B: 


Wu*ns = -W< 


• (0 X n^y (d • V Wg) • (13 + PgUn on and 


W, 


where (app. B, eq. (21)) 


(15) 


/ D A 

57 (^)i 


(16) 
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and (app. B, eq,(B.15) 



( 17 ) 


In addition, across wake surfaces the unsteady disturbance pressure must be continuous, i.e., 

[CpJ=0onW,, (18) 

when the unsteady disturbance pressure coefficient, i.e., 

is approximated as follows: 



-P U ^ 
2 0 ^0 



( 20 ) 


where eq.(B.29) has been introduced into eq.(19). Finally, the unsteady disturbance aerodynamic 
force is computed by integrating the pressure at the aerodynamic surfaces, viz.. 



when the pressure coefficient is given by eq. (20). 
Harmonic time dependence is imposed by substituting 


D = D* eioJt 


( 21 ) 


( 22 ) 


and 


0U " * 


eiwt 


(23) 


into the flow equation and boundary conditions. The result of that substitution into the 
flow equation shown as eq. (5) is as follows: 


V- (!>'■* = 2 ^ 2i£b (!>'* + 


(24) 
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where the freestream direction has been chosen to be sufficiently close to the direction of the x axis 
of the compressibility coordinate system that the components of 


,U„ 


• * - 0 X* 


have negligible order of magnitude. The result of substituting eq.(22) and eq.(23) into the 
boundary condition shown as eq. (15) is as follows t : 

w' * • rig = -Wg • (d* X tig^- (d* • vWg^ • rig + Pg Uj^i* on Sg and Wg 
where Wg (cf. eq. (10)) may contain a small component of freestream incidence, 

w'*= Pj,^v0'*-Mo^(io50' * +0(, D*), and un* = iw Uq D* ; 

while the wake boundary condition shown as eq. (12) becomes 

|[Cp*]=0onWg 


(25) 


( 26 ) 


(27) 


where 


n * zz 

P 1 


P U " 
2 « 0 


^Pg iw Uj, 0' * + Wg • v0'=*). 


A modified complex potential amplitude, , is introduced by letting 

0 / 'Ai 


0* ^ . 


(28) 


(29) 


Substituting eq. (29) into equations (24) - (28) leads to the following problem in which the 
flow equation assumes the form of Helmholtz’s equation, viz., 


<t>tz+ ^ ^ 

where 

n = o3 Mq/ 


t As shown in sec. 4.3, the method_for_c(^puting the steady component of flow is not sufficiently 
accurate for computing the term (D*‘ v Wg) • Ug. This term, therefore, is retained in the analysis but 
is not actually evaluated in the test cases of sec. 4. 
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and the boundary conditions appear as 


where 


* • ns =[- Wg • (r* X ng) -(d* • vWs) • % + un=^] on Sg and 


Ws (31) 


Further, eq. (28) becomes 
where 


Cp* = 0 on Wg 


(32) 

(33) 


p * = 


2 

a 


^ito / 0* + "^ Wg ‘ v0 


acjx 




(34) 


As will be seen in the following, a more useful form for the wake boundary condition follows by 
substituting eq. (34) into eq. (33). The desired form is given by 

ju* + — Wg • vp'* = 0 On Wg (35) 

where 


P* = lf0*]|. 


Letting 


a 1 - 

w — = — W„ • V 

P^ ^ 


(36) 


where £ is a coordinate line along a mean steady flow streamline on the wake surface and where 


Wc. 


-W 
P s 


is the magnitude of the nondimensional mass flux at the wake surface, one may write the following 
identity: 


(ico /p2 w,) M- + / “fAs ] 4^ (e 0“ A')). 
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( 37 ) 


Integrating along the steady streamline C from C = c, the boundary condition becomes 


M* (£, s) = ju* (c, s) e 


J d^/wg] 


where s is a wake surface coordinate which is orthogonal to £ . 

This expression evaluates the doublet strength at any point on a wake surface in terms of 
the value of the doublet strength at £ = c. As noted at the close of app. A, the doublet 
strength at the trailing edge of the lifting surface, from whence the wake emanates, is equal to 
the circulation on the lifting surface. This condition, when combined with that shown as 
eq. (37), allows the wake surface doublet distributions to be treated as known functions of 
the source and doublet distributions on the non-wake surfaces. This treatment of the wake 
doublet distribution is used in the following section. 


3.3 INTEGRAL SOLUTION TO UNSTEADY FLOW PROBLEMS 

As shown by app. C, Helmholtz’s theorem provides a solution to eq. (30), i.e., the flow equa- 
tion expressed in terms of the modified potential introduced by eq. (29). The solution provides 
the value of complex potential, 0*, at a point P interior to a volume V surrounded by the surface 
2, whereon unsteady sources and doublets are distributed. Letting Q represent a surface point 
of 2, the solution appears as follows: 


//f 

2 _ 


a* (Q) 0* (P, Q) + M* (Q) — - (0* (P, Q)) 


911;. 


ds' =0* (P) 


(38) 


where 




-ifiR 


R = |(x' - x')^ +j3^(y' - y)^ + j3^^z' - ^ , 


Q is a surface point with coordinates x ' , y ' , z ' , 
P is a spatial point with coordinates x , y , z , 
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K 


= nc • V 


and 


n' = |S2(n • t)t+ (n •?)) + (n • k) k . 


(39) 


(40) 


(41) 
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The vector n^, i$ the conormal vector of S defined such that, as shown by app. C, 



The function is an elementary solution to eq. (30) except at Q because * becomes singular 
as P approaches Q on the surface 2 . Multiplying this function by the factor employed to obtain 
the modified potential, viz.. 


e-i^Mg (x ' - x) , 

results in 

(p * = e''*^[^o(^ 

® 47tR ^ 

This equation appears as eq. (4.24) in ref. 6 where it is shown to represent the complex ampli- 
tude of the potential induced at point P due to a unit strength, harmonically fluctuating source 
located at the point Q. 

A solution to the boundary value problem posed in the preceding section, i.e., sec. 3.2, is con- 
structed by choosing the singularity distributions, a* (Q) and n* (Q), such that <l>* satisfies the 
flow boundary conditions. The boundary conditions are specified at the aerodynamic surface, 
Sg, the upper and lower surfaces of the wake, W^j and Wg, and, at the far-field surface, - 
the latter being an arbitrarily shaped surface located at a large distance from the aerodynamic 
surface. The original surface, 2, therefore is represented as 


2 = SsUWuUWgU2^ (43) 

where U denotes a union. 

The far-field boundary conditions require that the potential and its normal derivative vanish at 
These boundary conditions are satisfied by requiring a* (Q) and p* (Q) to vanish for points Q on 
2 qo ; as a consequence, the integral equation for 0* reduces to an integral over Sg, W^, and Wg. 


As shown in app. C, the source and doublet distribution strengths are related to the potential as follows; 


a* (Q) = 


d<p*/dn^ 


(Q) = [0*1 


(44) 
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where, for example, 


denotes the jump (or difference) in the limiting values of 0* as Q is approached from either side of 
the surface. Applying the second of eq. (44) to the wake surfaces leads to 

and 

«• (Qj) =<>♦ (Qi)-r (Qe)onWj. 

Letting the surface Wg approach the surface W^, in the limit 

0* (Qu) = </>* (Qc) 

SO that one may define 

(q;)-o* (Qj) 

as the jump in potential across the wake. Similarly, one may define 


[30*/ dn^] 3 (3^./ ^ - (30./3n,) ^ 

Qu Qg 


Further, noting that 



(45) 


the integrals over and Wg can be expressed as a single integral over the surface Wy, viz., 

// ([[90* / 9n^]]0 * + [[0*Ja0*/ 3nj,j ds . 
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The wake boundary conditions are shown by eq. (31). If eq. (31) is evaluated at either side 

of the wake surface and the evaluations are subtracted from one another, then we have the following 

result: 

|w**ns| = 0 on Wg . 


From eq. (32) and eq. (41) it follows that this expression may be written as follows: 


= 0 on Wg . (46) 

As a consequence of eq. (46), the first term of the wake integral integrand is seen to vanish. 

Using the above, the integral shown by eq. (38) becomes as follows: 


90* 

9nc 


0*(P) = 


//[‘ 


,9«c 


(P,Q) 


ds' + ^*(P,Q)ds' 


(47) 


Wc 


where the upper surface of the wake, , is identified as Wg , thereby distinguishing that side to be 
the side where n^ is positive. Also, as a consequence of eq. (37), eq. (47) becomes 


<p*(P) = 


//[” 


•(Q)**(P,Q) + m*(Q)^, 
9n . 


'I' (P.Q) 


ds' + 0^ (m* (c,s)) 


(48) 


where 


5 

<.w*(k*(c.s))= **(P,Q)ds' 

and ja* (c,s) denotes the doublet strength (or jump in potential, 10* 1) along the edge of the aero- 
dynamic surface from which the wake surface emanates. 


18 


3.4 AERODYNAMIC INFLUENCE COEFFICIENTS I 

The concepts involved in the panel method have been introduced in section 2.0, where it was indicated 
that the objective of the panel method is to reduce an integral equation and its associated boundary 
conditions to a determinate system of algebraic equations. The surfaces of integration are subdivided 
into panels and the unknown functions under the integral sign (viz., the singularity strength distribution 
functions) are approximated on each panel by polynomial functions of the surface coordinates % , t] . 
The coefficients of these polynomial functions are expressed in terms of a finite number of singularity 
parameters. The resulting expression is integrated over each panel to obtain influence coefficients 
which, when used with the boundary conditions, yield an algebraic system of equations containing 
the singularity parameters as the unknowns of the problem. 

The integral equation for the unsteady flow problem is shown by eq. (48) and the associated 
boundary conditions are given by eq. (31) and eq. (37). The panel method for reducing this 
problem to a system of algebraic equations is essentially that of ref. 1 . The paneling of the surfaces 
of integration and the panel approximations to the unknown surface functions are identical with 
those of ref. 1 . The kernel functions of the integral equation, however, are not identical with 
those of ref. 1 ; hence, different methods for evaluating the panel integrals are required. 

Regardless of this distinction there is a very close relationship with the panel method of ref. 1 . 

Because of this close relationship only the methods for evaluating the panel integrals are com- 
pletely developed in this report with the details of that development contained in app. D. 

This section presents only an overview of the panel method. 

The surface panels are defined by covering the actual surfaces of integration by a number of grids (cf., 
fig. 2) and these grids are called networks. The surface points at intersections of the grid lines (called 
grid points) are the quantities used to define the panels, cf. sec. 4.1 of ref. 1. The grid points of each 
network of panels are numbered by a double index system (M,N) which assumes that the grid points 
of a network are arranged in columns and rows even though the grid may not be rectangular, cf. sec. 1.1 
of ref. 2. If a surface of integration covered by a network has finite curvature, the panels are not 
actually segments of the curved surface; they lie close to the surface, but even their corner points need 
not coincide with the grid points, cf. fig. 6. The panels are flat or they are made up of several flat 
segments called subpanels; therefore, the panels are assembled to form a surface which is an approxi- 
mation to the actual surface of integration when that surface has curvature. 

Actually, depending on the distance from the panel to the point where its influence is being evaluated 
(i.e., the field point), one of three different surface of integration panel approximations is used. 

Letting D represent the panel diameter (viz., twice the distance from the panel center to its farthest 
corner) and letting Rq represent the distance from the panel center to the field point, the choice of 
panel approximation is as follows; 

(1) If Rq > 0.75D, then each panel is represented by the quadrilateral whose corners are 
formed by projecting four contiguous grid points to their average plane (ref. 1 , sec. 4.4). 

If the surface of integration has curvature, then the lines forming the edges of two adjacent 
panels need not coincide but need merely to intersect one another at their midpoints 
(cf. fig. 7). This approximate surface of integration, therefore, may be discontinuous 
at the panel edges. 
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/ 

Figure 6. — Panel Approximation to Surface of Integration 


/ 



Figure 7. — Approximate Surface Discontinuity 
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(2) If 0.45D < Rq < 0.75D, then the quadrilateral - viz., the four sided region formed by lines 

connecting four contiguous grid points - is represented by two plane triangular subpanels having 
vertices at the grid points and one edge formed by the diagonal, fig. 8. Since there are two 
diagonals, there are two possible surface approximations contained in this arrangement. Both 
approximations are used, in that, panel influences are computed using both surface approxima- 
tions and the quadrilateral panel influence is defined as the average influence. 



Figure 8. — Two Triangular Subpanels of a Quadrilateral Panel 


(3) If Rq< 0.45 D , then the quadrilateral is subdivided into eight triangular subpanels, fig. 9. 

The four interior subpanels lie in the average plane, while only the interior edge of an outer 
subpanel lies in the average panel plane. The exterior vertices of the outer subpanels are made to 
coincide with the grid points. This arrangement, like (2), above, causes the approximate surface 
of integration to be continuous at the panel edges. This continuity allows an approximate doublet 
distribution formulation having continuity at the panel edges, see sec. 4.2, ref. 1. 



Figure 9. — Eight Triangular Subpanels, Panel Center Local Coordinate System, 
and Nine Canonical Points 
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The panel approxima ting-functions for the source and doublet distributions are polynomial functions 
of local panel coordinate systems defined for each flat quadrilateral panel and each flat triangular 
subpanel. The coordinate lines of a local panel system are denoted as ^ , t? , f and the ^ , tj coor- 
dinate plane is chosen to coincide with the plane of a flat quadrilateral panel or a flat triangular 
subpanel. For case (1 ), above, the origin of the local coordinates is at the center of the quadrilateral. 
For case (2), above, the origin for each triangular subpanel is at the vertex opposite the side forming 
the diagonal of the quadrilateral. For case (3), the origin is at the center of the quadrilateral when the 
local coordinates are related to the inner four subpanels, fig. 9, and at the outer vertices when 
related to outer four subpanels, fig. 10. 



Figure 10. — Subpanel Approximation to Surface of Integration 


The polynomial approximating-functions for the source and doublet distributions are as follows: 
Source, 


a* (Q) = Uq* + + a^* 17 


(49) 


Doublet, 


M* (Q) - Mo* + M^* ^ ^ + 2^VV* ^ 


(50) 
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where Q is the surface point with coordinates t , t? . The values of the coefficients in these poly- 
nomials are found by fitting the polynomials to the values of the source and doublet distributions at 
the nine canonical points of the quadrilateral shown by fig. 9. The values at these nine points are 
denoted as 


and 


(e) 

»i* = a* (Qj) for i = 1, ... , 9 


(51) 


(e) 

Mi*=p*(Qi)for i- 1,... ,9. (52) 


For each panel and subpanel the coefficients are related to the canonical values by matrix expressions 
of the following form: 



Mq 






< 


Mr 


II * 


> 


k 


(e) 

(QQ)jik 



( 54 ) 


where the index j, in eq. (53) and eq. (54), ranges over the elements of the column matrices of 
polynomial coefficients, the index i ranges over the canonical values required to determine the 
polynomial coefficients, and the index k identifies the sub panel over which the approximating 
polynomial is defined. 
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The canonical values of the source and doublet distributions on a panel are determined from a 
least squares fit of the polynomials shown by eq. (49) and eq. (50) to the values of the distri- 
butions at the panel center and the centers of the neighboring panels. In the case of a panel 
source distribution the points of evaluation are those shown by fig. 1 1 , while in the case of a 
panel doublet distribution the points of evaluation are those shown by fig. 12. The distribution 
values at these points are referred to as the panel singularity parameters; and, for a panel source 
distribution, they are nine in number and denoted as 

(e) 

Uq,* for a= 1, ... , 9 . (55) 


For a panel doublet distribution the panel singularity parameters are twenty-one in number 
and are denoted as 

(e) 

Met* for « = 1 , ••• , 21. 



(56) 


Figure 1 1. - Source Evaluation Points at the Center of the Pane! and its Neighboring Panels 



Panel under consideration 


Doublet evaluation point 


Figure 12. — Doublet Evaluation Points at the Center of the Panel and Its Neighboring Panels 
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The least squares fit to these panel singularity parameters, when evaluated at the canonical points of 
the panel, yield the following matrix relations: 

j(e) ) r (e) 11(e) \ 

j aj*| = (ASTS)i„ (57) 

and 

f(e) if (e) ]f(e) | 

; (ASTD),^ (58) 

Taken together, eq. (49), eq. (50), eq. (53), eq. (54), eq. (57), eq. (58) describe the source and 
doublet distribution approximations on each subpanel of a quadrilateral, and this description 
is in terms of the panel source and doublet singularity parameters defined by eq. (55) and eq. (56). 

In the implementation of eq. (56) the points of evaluation appear to be 25 in number, because 
a double index system is used to enumerate the evaluation points. The double index system 
ranges over five columns and five rows of grid points in order to identify the twenty-one evaluation 
points. In other words, the evaluation points for each panel span five rows and five columns, 
and this leads to a single index ranging from one to twenty five. 

Having approximated the surfaces of integration by networks of panels and haying replaced the 
unknown functions describing the source and doublet distributions by approxirhating functions 
defined independently over each panel or subpanel, panel influence coefficients may be defined. 
Letting represent the surface area of the eth panel and letting Ug* and ju^* represent the 
approximating functions on the panel, the panel influences at the field point P are defined as follows : 

= (59) 

Sg 



( 60 ) 


( 61 ) 


(62) 
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where 


vi//* =-^' v//* ( 63 ) 

because, as shown by eq. (40),^* is a function of the coordinates of two points, P and Q, and the 
unprimed gradient is with respect to the coordinates of P while the primed gradient is with respect 
to the coordinates of Q. 


Appendix C contains the development whereby these integrals are expressed in the local panel 
coordinates and the results are shown as equations (C. 1 28) through (C. 1 36). The procedures 
for evaluating these integrals appear in app. D; they lead to a reduction of eqs. (59) through 
(63) to algebraic expressions evaluated at a finite number of points called control points, see 
sec. 3.4 and sec. 4, ref. 1 . These expressions are evaluated at control points located essentially 
at the singularity parameter evaluation points described above (see sec. D.l, app. D);if they 
are evaluated at the ith control point and combined with eqs. (49 through (58), they provide 
the following relations: 

9 3 (e) (e) (e) 

r; (DVDS*)jj (ASTS):„ 0 ^* 

= 1 j=l 


21 3 (e) (e) (e) 

(DVDD=^)ij(ASTD)j„p„* (65) 

ce=l j=l 

where the superscript (e) indicates the values of the quantities associated with the eth panel. 

The subscript i now indicates the ith control point. 



The indices a, appearing in eqs. (64) and (65), range over the singularity parameter evaluation 
points neighboring the eth panel. These equations and this index assignment scheme are set up 
exclusive of the arrangement of the panels when they are used to assemble an approximation to 
a surface of integration. A separate index assignment scheme is used to identify the singularity 
evaluation points on the actual surface (cf., secs. 4.1 and 4.2 of ref. 1 — particularly, the dis- 
cussion related to figs. 13 and 14). The terms “panel singularity parameters” and “global 
singularity parameters” are used. The values of the source and doublet distribution strengths 
appearing in eqs. (64) and (65) are referred to as panel singularity parameters; thus, a single 
isolated panel is under consideration and the panel singularity parameters at points on and 
neighboring that panel determine the doublet and/or source distribution on the panel. The 
values of the source and doublet distribution at the singularity evaluation points on the 
assembled surface are referred to as global singularity parameters. Incidence relations are used 
to relate the indices of the panel singularity parameters to those of the global singularity 
parameters. These relations appear as follows: 

(e) 

i = (llD)c^for a= 1 , ... , 21 (66) 


and 


(e) 

j = (IIS)q, for a = 1 . ... , 9 


(67) 
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where i is the global index corresponding to the ath doublet strength value for the eth panel and, 
similarly, j is the global index corresponding to the eth panel. These relations can be viewed as 
column matrices, one for each panel ; the ath entry of [(IID)a] is then the global index of a doublet 
evaluation point bn the assembled surface where the ath panel singularity point of the eth panel is 
to be located. 


The incidence relations shown by eqs. (66) and (67) are used to incorporate eqs. (64) and 
(65) into a single matrix expression, viz.. 



[(DSDFS.)^]{Xj*) 


( 68 ) 


where Xj* is the jth global singularity parameter. The coefficient matrix appearing in eq. (68) is 
called the aerodynamic influence coefficient matrix and this coefficient matrix contains the coefficients 
denoted as Ay* in sec. 2.0. The left hand member contains the complex amplitudes of the potential 
and the velocity (in terms of the modified potential, viz., eq. (29)) evaluated at the control points. 
However, as shown by app. C, sec. C.2, the potential and the velocity are discontinuous at the control 
points because the control points lie on surfaces containing the source and doublet distributions. The 
subscript A, appearing in eq. (68), denotes the fact that 
the average of the limiting values found by approaching 


the values computed by eq. (68) are 
the surface from above and below, viz., 

+ (69) 


and 


The discontinuities in potential and velocity at a surface of source and doublet distribution are as 
follows: 

For doublet distributions: 

(71) 

and for source distributions: 



I 
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( 73 ) 


thus, the potential and velocity at the upper side of the surface are found as follows: 


VM^ + rt^cr^, (74) 

the velocity at the upper side of the surface is found to be 

^u* = (■^)a + ^ ("c + VM*) ; 

and, at the lower side of the surface, similarly, 

0 ) 2 * = ( 75 ) 

■vg*=(^)A-Y(n^,a*+Vju*) (76) 

As a comment regarding eqs. {69) and (70), the average of the limiting values are not actually computed 
in the manner indicated by those equations. Although the integral describing the influence of a panel 
on a point located on the panel is singular, the singular behavior of the two integrals implied by the 
right hand numbers of eqs. (69) and (70) cancel one another. The result of this cancellation is a single, 
regular integral. This regular integral is evaluated to obtain the average influence, cf., ref. 1, sec. 3.7. 

3.5 UNSTEADY FLOW PROBLEM FORMULATION 

As a preliminary to this section, consider the network arrangement for describing the panels on the 
surface of a configuration (cf. sec. 4, ref. 1 and sec. 1.1, ref. 2). The panels are described in terms of 
networks consisting of rows and columns of panels. Each network covers a portion of the surface of a 
configuration without overlapping the surface covered by another network; thus, every panel on the 
surface of a configuration is a member of one and only one network. For example, the right half of 
a symmetric wing surface and its wake could be covered by four networks consisting of the upper 
wing surface, the lower wing surface, a surface closing the tip, and the wake surface. 

The panels of a single network are all of one type; that is, every panel of a particular network is either 
a source panel, a doublet panel, or a combined source-doublet panel. If one portion of a configuration 
surface is to have only sources distributed on it while a second portion is to have only doublets dis- 
tributed on it, then these two portions of surface must be covered by different networks: one a source 
network, the other a doublet network. 

As already noted, a network is a two dimensional array of panels, fig. 13. The grid points, which form 
the comers of the panels of a network, are enumerated by the indices (M,N) where M is the grid 
point row number, while N is the grid point column number. Although it would appear from this 
that a network always has four sides, one side of a network can be collapsed into a single point so 
that the collapsed side has zero length. When the edge (or side) of a network is collapsed, all panels 
along that edge must have a collapsed edge such that they are three sided and have one corner point 


0u* = (<^*)a 


Noting that 


v^* =Vju* -I- 


'30*' 

9n„ 


(30* ' 
i9n^ , 
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located at the point representing the collapsed edge of the network. This requirement is consistent 
with the requirement that every row and column of grid points of a network must have the same number 
of members. If the collapsed edge of a network corresponds with the column direction of indexing 
grid points and the network has M rows of grid points, then the point representing the collapsed edge 
is counted M times in the indexing scheme. Also, if a network edge has finite length, then the edges 
of every panel lying on that network edge must have a finite length even though the length of each 
panel edge can be different from the others. 

Recall from the preceding that the unit vector normal to the surface appears in the boundary conditions 
of the flow problem; thus, in order to formulate a flow problem one must have a facility for specifying 
which side of a surface is to be that of positive normal (i.e., the upper side of the surface). This facility 
is provided for each network by the grid point numbering system. Let M and N denote vectors, 
respectively, in the grid point row and column directions (fig. 13). The side of the surface in the 
direction of the cross product (N x M) is called the upper side (or positive side) of the surface and 
this is the side of positive unit normal vector for the portion of surface covered by the network. The 
opposite side is called the lower (or negative) side. 

Network surfaces do not overlap one another and the edges of networks are either free edges of surfaces 
or edges which abut other network edges. Since surfaces can contain narrow gaps, a convention is used 
to distinguish between a physical gap and an abutment of two network edges on a continuous surface. 
The convention is as follows: The edge of a network is a continuous curve consisting of straight line 
segments. If the two curves representing the edges of two networks are identical (within the accuracy 
of the computations), then the two network edges abut one another. Along any straight line segment 
of the network edge curves, the number and location of the grid points can differ between the abutting 
networks (except at the line segment ends) but these intermediate grid points must lie precisely on 
the line segment, fig. 14. If these conditions are not satisfied, then the two networks do not abut one 
another and a surface gap (or discontinuity) is introduced. 

The control points of typical networks are shown by fig. 15. A source network has panel center 
control points only ; while a doublet network has both panel center and network edge control 
points.f Boundary conditions are applied at these points specifying: (1) the value of a flow 
parameter obtained as a value of an aerodynamic influence coefficient (e.g., the normal compo- 
nent of the unsteady perturbation mass flux vector at the upper side of the network) or (2) the 
value of a source or doublet distribution. The network edge control points play the role of 
outer evaluation points, see app. C of ref. 1 ; they are used to impose continuity on the doublet 
distribution between abutting networks of doublet panels, see sec. 6.0, ref. 1 . 

The aerodynamic surface boundary condition which appears as eq. (31) and which is derived in app. B 
is expressed as follows: 

w* • Og = B* on Sg 


where 


B= 


= [- Ws • (« ‘x %) - (d’ 


VWj) 


I 

Hs+P, 


s n 




(31) 


t Wake surfaces are represented by doublet networks; and, in the case of a wake-doublet network, 
control points occur only along one edge, cf., the discussion concerning figures 16 and 17. 
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Source network 


Doublet network 


Figure 15. — Network Control Points 


This boundary condition must be satisfied at each of the aerodynamic surface network control points. 
Referring to eqs. (32) and (41), the left hand number of eq. (31) can also be expressed in 
terms of the mean steady surface conormal vector such that 


n„ • V0* = B* on S„ ; 
0 


(77) 


hence, assuming that this boundary condition is to be satisfied at the upper surface of a network, it 
is evaluated, using eq. (72), as follows: 


1 ^ 1 ^ \ 

Vy^* +-n^a* +-VJU*) on . 


(78) 


Eq. (78) represents a special form of the following general boundary condition: 

C„(w„*-a)+f„-V + D„V+Ce(w«*-8)+T8-V + Dj0j* = B«. (79) 

If the coefficients of the flow parameters are given the following values: ' 

C„=l, T„ = fj = C, = D„=D(,=0, 

then eq. (79) reduces directly to eq. (78). 

Consider, as an example, the Morino-type boundary conditions (ref. 2, sec. 1.31) for an aerodynamic 
body having finite thickness, i.e., a body whose surface encloses a volume of space. These boundary 
conditions are such that the unsteady disturbance potential is made to vanish everywhere in the enclosed 
volume. The networks on the surface are chosen as combined source-doublet networks and their 
grid point indices are chosen such that the lower side of the network surfaces are at the interior of 
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the body. If (j>^* is set to zero at the doublet control points, then 0* will vanish at the lower side of 
the surface; and, therefore, <j>* will vanish everywhere throughout the enclosed volume. Since, as shown 
by eq. (71 ), the doublet strength is equal to the jump in potential at the surface, i.e., 

/x* = 0u*-0e*, 

the doublet strength is equal to the potential at the upper side of the surface, viz., 

M* = 0u* ’ 

because 

0g* = 0 . 

• Also, from cq. (72) 

'a0*V 

but, because 

£ 

the source strength is equal to the normal component of pertqrbation mass flux at the upper side of the 
surface, viz.. 


or 

o*= (^n • w*)y 

where eqs. (32) and (41 ) have been used to introduce the mass flux vector. 

The Morino-type boundary conditions are implemented, using the general form of the boundary 
condition shown as eq. (79), by the following choice of coefficients; 

At doublet control points 





set Dg = 1 and all other coefficients set to zero. 

At source control points 

set Cy = 1 and set B* = (right hand side of eq. (31) ) 
and all other coefficients set to zero. 
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The evaluation of the right hand side of eq. (31) requires the values of the mass flux vector, the mass 
flux vector gradient, and the mass density of the fluid in the steady mean component of flow. These 
quantities must be obtained from a separate solution to the steady mean flow problem, cf. sec. 5.2. 


In the case of thin wing theory (ref. 7, Chapter 13) the steady mean component of flow is simply 
the undistribed freestream. In that case the aerodynamic surface boundary condition, viz., eq. (31), 
reduces to the following: 


w^ 


A 


A 

-1 


(« • X "s) 


+ u„ V U„ 


g-icoxMo^/ j3^ 


on S, 


(80) 


where is a cylindrical surface with generator parallel with the x axis. Only doublet networks are 
required on Sg and the boundary condition is expressed by setting Cg = 0.5, C^ = 0.5, and B* 
equal to the right hand member of eq. (80) evaluated at the doublet control points. All other 
coefficients of the boundary condition are set to zero. 

As indicated by eq. (48), the unsteady wake surface boundary condition is not an independent 
boundary condition, in that the doublet distribution on the wake is determined on its entire 
surface by its value along its upstream edge. The wake is a surface emanating from the edge of 
a lifting body, fig. 16, and the Kutta condition requires that the discontinuity in potential across 
the wake surface and across the lifting surface be continuous at the lifting surface trailing edge. 

The wake is represented by doublet networks; and, since the doublet strength is equal to the 
jump in potential, i.e.. 


M* (Q) = 0u* (Q) - h* (Q) Q o" Wg , 

the wake doublet strength along the network edges which abut a lifting surface must be equal to the 
jump in potential across the lifting surface along the line of abutment. If the lifting body is 
represented by thin wing theory, then the lifting body is an infinitesimally thin lifting surface 
represented by doublet networks. In this case, letting Qg represent a point on the lifting surface and 
letting represent a point on the wake, if Qg and approach the same point on the abutment line 
between the lifting surface and the wake, then in the limit we require that 


M* (Qs)"/^* (^w) > 


(81) 



Figure 16. — Surfaces of a Lifting Body and its Wake. 
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and this expression represents the Kutta condition requirement. If the lifting surface is 
thick and is represented by Morino-type boundary conditions such that 




and 


M* (Qs^) = (Qs^) 

where and Qg represent points on the two surfaces shown by fig. 16, then in the limit 
as the three points Qg , Qg^ and all approach the same point on the line of abutment 
with the wake we require that 

M* (q/) -M* (Qs^) (Qw) ; 

and this expression, in addition to eq. (37), represents the Kutta condition requirement.t 


(82) 


As noted above, the doublet strength distribution on a wake network is determined by the conditions 
along the upstream edge of the network. This determination is described by eq. (37), viz.. 


jU* (fi , s) = JU* (c , s) e 


(ico /<3^) r (C , s) 


on W., 


(37) 


where 


C 

T (£ , s) = y* 

c 


Wg (^ , c) 


Po^o 

is the magnitude of the mean steady mass flux vector divided by the freestream mass flux. The parameter 
8 is a coordinate line along the mean steady flow streamline at the wake surface, s is a coordinate 
line orthogonal to 8 , and c is the value of 8 at the upstream edge of the wake network, see fig. 1 7. 

Eq. (37) can be seen to specify the doublet distribution everywhere on a wake network in terms of 
the doublets’ strength along the upstream edge. 

If, as in the case of the usual linear theory, the mean steady flow is simply the undisturbed 
freestream, then eq. (37) reduces to 



tx* (x , s) = p* (c , s)e (83) 


+ Asa result of numerical approximations involved in the evaluation of eq. (82), a modified 
form of the Morino-type boundary conditions is required for approximating the Kutta 
condition, cf., sec. 4.4. 
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Figure 17. — Surfaces of a Lifting Body and Its Wake 


the X coordinate. The results computed in the following section have all been computed using 
this doublet distribution on the wake surface networks. The theory and its implementation as 
a computer program provide for the wake model shown as eq. (37), but the steady flow analysis 
capability for generating the parameters required to evaluate eq. (37) has not yet been developed. 

The evaluation of eq. (37) requires a solution to a nonlinear steady flow problem. The steady flow 
method of solution provided by ref. 1 is based on a wake which is a cylindrical surface with gen- 
erator parallel with the undisturbed freestream direction. This wake representation provides a 
first order approximation to a solution satisfying the exact wake boundary conditions developed 
in app. A, viz., 


W'OonW (A.35) 

and 

Vj^ • n = 0 on W ; (A.48) 

however, it does not provide the data required to evaluate eq. (37). 

In the exact theory eq. (A. 35) requires pressure to be continuous across the wake and eq. (A.48) 
requires the mean velocity at the wake surface to be tangent to the wake surface. A plausible 
formulation of these boundary conditions, which is consistent with the approximations involved 
in eq. (37), consists of replacing the pressure in eq. (A. 35) with the second order approximation 
shown by eq. (B.8) and of replacing the flow velocity appearing in eq. (A.48) with the approximate 
mass flux vector shown by eq. (B.4). An iterative solution is required in which the wake surface 
location is adjusted until these approximate forms of the wake boundary conditions are satisfied. 
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Presuming success for the solution of the problem, the solution would provide the spatial coordi- 
nates of the wake surface (viz,, £, s) and the mean of the steady mass flux vector along the wake 
surface (viz., Wg). These quantities then provide a means for evaluating the unsteady flow wake 
boundary condition shown as eq. (37). 

Regardless of how the wake surface is approximated, the doublet strength varies along a streamline 
coordinate deviating only slightly from a harmonic function of the coordinate distance. Recalling 
that the panel approximation to the doublet strength distribution is a quadratic function of the local 
panel coordinates, cf. eq. (50), it can be seen that the wake paneling must be sufficiently dense 
in the streamwise direction that the panel-wide quadratic provides an accurate least squares 
approximation to the harmonic variation. A panel density of sixteen panels per wave length 
of harmonic doublet variation has been used in the example cases described in sec. 4. 
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4.0 DISCUSSION OF RESULTS 


This section presents computed results which demonstrate the validity of the approximations 
contained in the unsteady flow panel method, primarily, the approximation to the unsteady kernel 
function described in app. D. This emphasis on the validity of the kernal function approximation 
stems from the fact that the method is an outgrowth of the steady flow panel method of ref. 1 . 

The approximation to the unsteady kernel function is the primary feature leading to this outgrowth; 
hence, the demonstration of its validity is a central consideration in establishing the validity of 
the unsteady flow panel method. 

The demonstration cases consist of configurations of four different types: 

(1) thin, planar lifting surfaces; 

(2) thin, nonplanar lifting surfaces; 

(3) bodies of finite thickness; and 

(4) combinations of (.2) and (3). 

Demonstration cases of type (1) are the primary basis for validating the unsteady kernel function 
approximation. The reason for this is the availability of an accurate, alternate method as a basis 
of comparison, viz., that of ref. 8. A T-tail is typical of a type (2) configuration and three different 
T-tails are evaluated and compared with alternate methods of evaluation. A type (3) configuration 
is represented by a wing with sources and doublets distributed on its actual surfaces in their mean 
steady locations. The lift distribution on this wing as a thick body is compared with the lift 
distribution on the same wing treated as a thin body, type (1). A type (4) configuration is repre- 
sented by a twin engined transport configuration. 


4. 1 APPLICATION OF PANEL METHOD TO PLANAR, THIN WINGS 

As noted in the preceding, in the case of thin lifting surfaces, when the boundary conditions are 
completely linearized, the steady mean component of flow consists of the undisturbed freestream. 
The aerodynamic surface and wake boundary conditions shown by eqs. (31) and (35), reduce 
to those shown by eqs. (80) and (83); and the panel method reduces to a method for evaluating 
the flow governed by thin wing theory, e.g., chapter 13 of ref. 7. The problem is formulated 
entirely in terms of doublet networks on the lifting surfaces and on the surfaces of their wakes. 
For the case of a single, planar lifting surface the acceleration doublet method of ref. 8 provides 
a well validated, independent method for evaluating thin wing theory; hence, it provides a 
basis for testing the validity and accuracy of the panel method doublet network for in-plane 
computations. 

Figure 1 8 shows the planform of the right hand wing of a test case wing, which is symmetric 
about the x axis. The wing is oscillating in pitch about the midpoint of the root chord with a 
reduced frequency of 0.3577 in a uniform flow with a Mach number of 0.6. Figures 19 through 
22 show the real and imaginary parts of the lifting pressure coefficient chordwise distribution at 
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Lifting pressure coefficient 



Figure 19. — Thin Wing Unsteady Lifting Pressure Comparison at 18. 1% Semi-span 
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Lifting pressure coefficient, ||C 




Panel method 
O real 

□ imaginary 
Reference 8 
O real 

<!> imaginary 


Fraction of chord, X/C 


Figure 20. — Thin Wing Unsteady Lifting Pressure Comparisons at 5 1.2% Semi-span 
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Lifting pressure coefficient, [[C 







four spanwise locations (viz., at 18.1, 5 1.2, 81.7, and 97.7 percent of semi span from the root). 

The plotted values shown on the figure as circles and squares are, respectively, the real and 
imaginary lifting pressure coefficients computed by the panel method using 100 panels on the 
right hand wing. The panel spacing is such that there are ten panels in both the spanwise and 
chordwise directions; letting b = 4.05 denote the span, the chordwise edges of the panels are located 
at the following points along the y coordinate direction: 


where 


b /S7T\ 

y(N)=j sm(-) 

(84) 

s = (N- 1)/(NMAX- 1) 

(85) 


and NMAX = 1 1 while N = 1 ,—, NMAX. The x coordinates of spanwise panel edges are spaced 
along each chord lines as follows: 

I (M)/c ~ ^ 

where c is the local chord length, § is the distance from the trailing edge, and 

s = (M- 1)/(MMAX-1) (87) 

where MMAX = 1 1 and M = 1,—, MMAX. Fig. 23 shows the resulting panel spacing. 



43 


Figures 19 tlirough 22 show the panel method predicted pressure distribution compared with 
that predicted by the method of ref. 8. This method is based on a solution to the unsteady 
flow problem in terms of the acceleration doublet kernel function and distribution functions 
defined over the entire wing planform; it is a sufficiently well validated method that the close 
agreement shown between it and the panel method validates the panel method for application 
to planar, lifting surfaces. 

Table 1 shows the effect of panel density on the values of the predicted unsteady pressures acting 
on the test wing. The effect is evaluated on the basis of the errors in the unsteady lift, pitching 
moment and rolling moment. This error evaluation is expressed as the relative error in the magni- 
tudes of the lift and the moments, and the absolute error in their phase angles. When the wing 
has 375 panels, viz.. 25 chordwise panel spaces and 15 spanwise panel spaces, the panel method 
yields the following: 

Lift Force: 

ICl*I = ' ^ ^ = 3.6235 

9Sr|hP 


where 


IL 


I is the magnitude of the complex lift force on the right hand wing, Sj^pp - 2.5465 


square semi root chord lengths, and q is the dynamic pressure. The lift phase angle is 

(Pp= 17.099° 


Rolling Moment: 


* 

1C, I 


lg*l 

9 ^RFF ^REF 


1.5672 


where IC I is the magnitude of the complex rolling moment on the right hand wing about the 
x-axis and Bj^pp = 1.9463 semi root chord lengths. The rolling moment phase angle is 

14.852° 

Pitching Moment: 


9 Crpp Srp^ 

where Im'^ I is the magnitude of the complex pitching moment about the pitch axis at x = 1 .0 
and Cjppp = 1.308 semi root chord lengths. 

The amplitude and phase errors are evaluated relative to tlie values of these quantities computed 
for the wing with 375 panels. 
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Table 1. — Effect of Panel Density on Thin Wing Forces 


Panel density 

Chordwise Spanwise 



Lift coefficient 

Roll coefficient 

Relative 

Absol. 

Relative 

Absol. 

amplitude 

phase 

amplitude 

phase 

error 

error 

error 

error 

(percent) 

(deg) 

(percent) 

(deg) 

1.2882 

3.7643 

-4.3078 

3.5575 

.0672 

2.3262 

-3.3658 

1.3043 

- .5647 

1.7591 

-2.9070 i 

.5999 

- .7855 

1.2611 

-2.7774 

.2514 

- .3229 

.9876 

-2.7966 

.2817 

.0096 

1.7587 

-1.2023 

1.0720 

- .5304 I 

1.1990 

-1.1313 

.4498 

- -8324 

.7316 

-1.1493 

.0848 

- .5755 

.5154 

-1.0051 

.0396 

’ .9754 

1.0403 

-1.5814 

.5754 

-1.3384 

.5922 

-1.7217 

.2047 

-1.1723 

.4087 

-1.5706 

.1310 

- .9496 

.3714 

-1.4320 

.1594 

-1.1424 

-1.0397 


-1.4593 

-1.3294 

.2290 

.1242 

- .8520 

.2100 

-1.1682 

.1299 

■ .7375 

.2077 

-1.0765 

.1425 

- .5940 

.0770 

- .7108 

.0478 

- .4169 

.0572 

- .5361 

.0401 

- .3069 

.0576 

- .4276 

.0463 

- .1727 

- .0043 

- .1689 

- .0084 

- .0661 

- .0049 

- ,0604 

- .0022 

0.0000 

0.0000 

0.0000 

0.0000 


Pitch coefficient 
Relative Absol. 

amplitude phase 

error error 

(percent) (deg) 


-1.4574 1.3216 

-1.0146 - .0571 

-3.2965 - .8675 

-3.0980 - .9688 

- .0434 .0365 

-6.9601 -4.7459 

-6.1607 -3.1378 

-5.4025 -2.4572 

-2.5278 -1.1698 

-7.4157 -3.8772 

-6.1414 -2.6490 

-3.0597 -1.0612 

-1.0267 - .1956 

-6.3652 -2.8913 

-3.1733 -1.1056 

-1.2578 - .2654 

- .3041 .1269 

-3.2516 -1.3135 

-1.3554 

- .4553 
-1.3835 

- .4829 

0.0000 













Figure 24 shows a second test wing. This is a swept wing having no taper, a twenty-five degree 
angle of sweep, and partial span trailing edge flaps having chord lengths equal to thirty per cent 
of the total wing chord. This wing was tested for various combinations of flap and wing 
oscillatory motions and the results of those tests are reported by ref. 9. The real and imaginary 
parts of the unsteady pressure coefficient were measured along the chord lines shown as dashed 
lines on fig. 24. The measured pressure coefficient distributions are compared with distributions 
computed by the panel method in fig. 25 through 27. The surface paneling arrangement is 
shown in fig. 28. 

The flow velocity is less than 100 feet per second; hence, the flow was taken to be incompressible. 
There are five different cases of flap and wing oscillation, which are as follows: 

( 1 ) Both flaps oscillating with the same phase and amplitude (viz., 0.82 degrees) 

(2) Outer flap oscillating (amplitude: 0.66 degrees), inner flap and wing stationary. 



Figure 24. — Test Wing with Partial Span Flaps 
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Lifting pressure coefficient 
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Fraction of chord, X/C 


a. REAL PART 

Figure 25, — Chordwise Distribution of Lifting Pressure Coefficient, Case 


Experimental values 
0 -0.25m inboard 
□ -0.05m outboard 
0 -0.28m outboard 




0.0 .2 .4 .6 .8 1.0 

Fraction of chord, X/C 

b. IMAGINARY PART 

Figure 25. — Chordwise Distribution of Lifting Pressure Coefficient, Case 1. 


Experimental values 
O -0.25m inboard 

□ -0,05m outboard 
<!> -0.28m outboard 





Experimental values 
O -0.25m Inboard 

□ -6.05m outboard 
O -0.28m outboard 


Fraction of chord, X/C 


a. REAL PART 


Figure 26. — Chordwise Distribution of Lifting Pressure Coefficient, Case 2. 





Lifting coefficient 


.012 



b. IMAGINARY PART 

Figure 26. — Chordwise Distribution of Lifting Pressure Coefficient, Case 2. 
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(3) Inner flap oscillating (amplitude: 0.67 degrees), outer flap oscillating (amplitude: 

0.65 degrees) in anti-phase with inner flap, wing stationary. 

In each case the reduced frequency of the oscillatory motion is 0.372, 

The comparisons shown by fig. 25 through 27 demonstrate the capability of the panel method 
to evaluate the lifting pressure distributions induced by oscillating trailing edge flaps. Table 2 
further demonstrates that capability; it shows the variation in the amplitude of the control 
surface hinge moments for several panel densities. This variation is shown as the error in hinge 
moment amplitude relative to that computed for the most densely paneled case. This table 
also shows the absolute error in predicted phase angle, again, regarding the phase angle of the most 
densely paneled case as the point of reference. 



Experimental values 
O -0.25m inboard 
A -0.05m outboard 
IP -0.28m outboard 


(a.) REAL PART 


Figure 27. — Chordwise Distribution of Lifting Pressure Coefficient, Case 3 




Lifting coefficient 



(b,) IMAGINARY PART 

Figure 27. - Chord wise Distribution of Lifting Pressure Coefficient, Case 3 


Experimental values 
O -0.25m inboard 
□ -0.05m outboard 
<J> -0.28m outboard 
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Table 2. - Swept, Thin Wing Hinged Flap Hinge Moment Amplitude and Phase Errors 


Panel 

Chordwise 

M 

density 

Spanwise 

N 

Mode shape 
number 

Inboard flap 

Relative Absol. 

amplitude phase 

error error 

(percent) (deg) 

Outboard flap 

Relative Absol. 

amplitude . phase 

error error 

(percent) (deg) 

3 

3 

1 

8.3573 

- 1.0627 

7.0874 

- 2.3467 

4 

4 


6.8592 

- .5331 

3.6825 

- 1.4808 

5 

5 


3.941 1 

- .3346 

1.9150 

- 1.0299 

7 

7 


1 .6780 

- .1575 

.4573 

- .3566 

10 

7 


.0910 

- .0180 

- ,6030 

. 148 r 

10 

10 


0.0000 

0.0000 

0.0000 

0.0000 

3 

3 

2 

21.9326 

-17.4810 

6.7727 

- 2.7354 

4 

4 


8.9168 

- 10.3497 

3.6726 

- 1.7837 

5 

5 


4.0357 

- 6.1419 

1.7157 

- 1.1981 

7 

7 


1.4196 

- 2.0295 

.3400 

- .4614 

10 

7 


.3209 

- 1.1490 

- .7575 

.1303 

10 

10 


0.0000 

0.0000 

0.0000 

0.0000 

3 

3 

3 

10.6967 

- 1.6447 

6.2658 

- 3.2233 

4 

4 


7.9253 

- .9786 

3.5864 

- 2.1791 

5 

5 


4.5119 

- .6171 

1.4140 

- 1.4031 

7 

7 


1.8492 

- .2542 

.1613 

- .5891 

10 

7 


.2118 

- .0640 

- .9636 

.1196 

10 

10 


0.0000 

0.0000 

0.0000 

0.0000 




Figure 28. — Paneling of Test Wing with Split Trailing Edge Flap 


4.2 APPLICATION OF PANEL METHOD TO NON PLANAR , 

THIN LIFTING SURFACES 

The panel method is validated for this type of configuration by computing the unsteady aero- 
dynamic forces on three different T-tail configurations undergoing oscillatory motions. These 
are T- tails which have been evaluated by the alternative methods described by ref. 10 through 13. 

The first T-tail evaluated is shown by fig. 29. This T-tail is tested oscillating in yaw about an 
axis along the mid chord line of the vertical surface and the test results are contained in ref. 10. 
The tail oscillates with a reduced frequency of 0.1 05 in a uniform freestream having a Mach 
number of 0.376. Table 3 shows the values of the amplitude and phase angle of the following 
generalized aerodynamic force coefficients from three sources: 

Cy side force coefficient due to yaw 

C'g rolling moment coefficient of the horizontal surface abouts its center 

line due to yaw 

rolling moment coefficient of the T-tail about the root of the vertical 
due to yaw 

Cn^‘ yawing moment coefficient due to yaw 
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Figure 29. — Clevenson-Leadbetter T-tail 


One set of values appearing in table 3 is that from the experiment of ref. 9, the second set 
was computed by the method of ref. 1 1 , and the third set was computed by the panel method 
using the paneling shown by fig. 30. 
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Figure 30. — Paneling of Clevenson-Leadbetter T-Taii 


Figure 3 1 shows the T-tail of ref. 1 1 . The side force, yawing moment, and rolling moment coeffi- 
cients for this tail are evaluated for five different modes of motion. The five modes of motion are 
described in terms of mode shapes expressed in terms of the coordinate system shown on figure 3 1 
For points on the horizontal they are as follows: 

Zi=0 
Z 2 = 0 

Z3 = y 

24 = 0 

Z 5 -2 (1.5041)y 
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For points on the vertical they are as follows: 
yi = l 
y2 = x 

y3 = (1.5041 -Z) 
y4 = X(1.5041 -Z) 
y5 = (1.5041 -Z)2 

The oscillatory motions all have a reduced frequency of 0.5 based on the root chord length of the 
tailplane, and the flow Mach number is 0.8. 


D 



Reference 11 



Figure 31. — Plan forms of Horizontal and Vertical Fin of Davies T-taH 


58 


Table 4 lists the values of the real and imaginary parts of the three generalized aerodynamic 
force coefficients corresponding to the first, second, and third mode shapes, viz., Cy, Cj^, and Cg 
with the yawing moment, n, taken about the z axis shown on fig. 3 1 and the rolling moment, 9 . , 
taken about the root chord line of the fin. Complex values for these three coefficients relative 
to each of the above five modes of motion are shown by table 4 and the values shown are those 
presented by ref. 1 1 along with the corresponding values computed by the panel rnethod. The 
panel arrangement is essentially that shown by fig. 30 with 1 1 panel spaces in the spanwise 
direction and 13 panel spaces in the chordwise direction of each of the three lifting surfaces making 
up the T-tail. Table 5 shows the effect of varying the panel density on the predicted aerodynamics of 
this T-tail. Table 5 shows the relative error in the amplitudes and the absolute error in the phase 
angles for four different panel densities, treating the most densely paneled case as the basis. The res- 
ults of ref. 1 1 are included also, again, using the most densely paneled as a base for comparison. 


Table 4. — Comparison of Davies T-tail Generalized Force Coefficients 


Coefficients 

Mode 1 

Mode 2 

Mode 3 

Mode 4 

Mode 5 

Real 

Imag 

Real 

Imag 

Real 

Imag 

Real 

Imag 

Real 

Imag 


.2164 

.2236 

-1.0246 

-1.0490 

-2.0977 

-2.2056 

-.7804 

-.8461 

.1589 

.1431 

-.9654 

-.9327 

-1.700 

-1.8105 

-.6946 

-.7158 

.1529 

.1058 

-1.2249 

-1.0671 

S"' . Ref. 11 

Vawmg 

moment 

-.0129 

-.0327 

.3185 

.3049 

.6847 

.6776 

-.1498 

-.2869 

.0261 

.0673 

.1489 

.1274 

.3520 

.3361 

-.1612 

-.2926 

.0529 

.0937 

.0703 

.0317 

Ref. 11 

Rolling p^, 

moment 

.0859 

.0872 

-1.0942 

-1.0582 

-2.2202 

-2.2048 

-.6104 

-.6527 

.1759 

.0761 

-1.6627 

-1.3939 

-2.0574 

-2.0616 

-.7025 

-.6969 

.3505 

.0641 

-3.0422 

-2.2815 


* P.M. — Panel method computed result 


Figure 32 shows the T-tail which is evaluated by ref. 12 for oscillatory sideslip, yaw, and roll. 
LetL* denote the following generalized complex aerodynamic force coefficients: 


47rqb*^^M 


H^ds 


where for ju = 1, 2, 3 describes the T-tail surface displacement for unit yaw, sideslip, and 
roll while H*, for p = 1 , 2, 3 describes the complex amplitude of the surface distribution of 
lifting pressure induced by unit yaw, sideslip, and roll. Table 6 compares the values of these 
generalized force coefficients computed by the method of ref. 1 2 with the values computed 
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Table 5. — Davies T-tail Generalized Force Coefficient Amplitude and Phase Errors 


Panel density 


Mode shape 

Side force coefficient 

Roll moment coefficient 

Yaw moment coefficient 

Chordwise Spanwise 

Number 

Relative 

amplitude 

Absol. 

phase 

Relative 

amplitude 

Absol. 

phase 

Relative 

amplitude 

Absol. 

phase 




error 

error 

error 

error 

error 

error 

H 

N 


(percent) 

(deg) 

(percent) 

(deg) 

(percent) 

(deg) 

4 

4 

1 

• 7113 

1.7037 

6.7321 

2.0443 

-13.7859 

2.6247 



2 

2.1025 

1-2369 

3.1898 

1.6755 

-8.2479 

-340742 



3 

6.0527 

-.4840 

19.5132 

1.6025 

-17.9444 

-17.0173 



4 

3.3682 

.6625 

9.1164 

1.2442 

-3.6290 

-7.8951 



5 

13.6220 

-2.9102 

29.5909 

1.7166 

25.8300 

-37.0972 

10 

4 

1 

-2.3728 

-.0500 

2.3127 

.0467 

-4.2382 

• 9793 



2 

-2.2320 

-.2738 

2.5290 

-.3496 

-3.4322 

-.9620 



3 

2.9433 

-1.8463 

17.0718 

-.2244 

-11.4497 

-5.3214 



4 

-2.5776 

-.5884 

1-1952 

-.5656 

-1.6913 

-.7936 



5 

11.5950 

-4.0580 

29.6109 

-.0907 

-.2585 

-31.5162 

10 

9 

1 

.3192 

-.0389 

1.7467 

-.0757 

-1.8769 

.3091 



2 

.3685 

-.0522 

1.9557 

-.1367 

-1.3871 

-•4019 



3 

1.3557 

-.2629 

4.4635 

.0119 

-3.1876 

-1.6598 



4 

.7066 

-.0961 

1.9844 

-.1043 

-.8659 

-.9761 



5 

2.6007 

-.8238 

6-5778 

.1386 

-.6875 

-6.1689 

13 

11 

1 

O.OQOO 

0.0000 

0.0000 

0.0000 

8.0000 

0.0000 



2 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 



3 

0.0000 

0.0000 

0.0000 

0.0000 

0.0008 

0.0080 



4 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 



5 

0.0000 

0.0000 

Q.OOOO 

0.0008 

0.0000 

0.0000 

(Reference 11) 


1 

-2.3647 

-.1070 

3.3698 

-.2220 

3.9500 

-3.8021 



2 

-5.2558 

-.5810 

• 1391 

-1.1181 

-4.7433 

10.6073 



3 

3.6854 

• 6242 

19.7713 

2.9140 

4.9185 

17.9035 



4 

-5.6725 

• 6525 

-.1000 

.1753 

-13.1200 

16.4364 



5 

15.1142 

1.4530 

34.1712 

4.9629 

-11.0567 

34.3475 


b 



Reference 12 


Figure 32. — Plant arms of Horizontal Tail and Vertical Fin of Stark T-taU 


using the panel method. The panel method values are based on a computation using the panel 
arrangement shown by fig. 33. The comparisons are shown in table 6 for two different values 
of reduced frequency of motion, viz., 0.2 and 0.3. These reduced frequency values are based 
on the definition of ref. 1 2 using one third of the semi-span of the horizontal tail as the 
characteristic length. Also included in table 6 are comparison values computed using the 
doublet lattice method; those values are taken from ref. 13. 
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Figure 33. — Paneling of Stark T-tan 


Table 6. — Comparison of Stark T-taH Generalized Force Coefficients 


Coefficients 


Ref. 12 
P.M. 
Ref. 13 

Ref. 12 
P.M. 
Ref. 13 

Ref. 12 
P.M. 
Ref. 13 


Ref. 12 
P.M. 

Ref. 12 
P.M. 

Ref. 12 
P.M. 


Yawing 

Real Imag 

Sideslipping 
Real Imag 

Rolling 
Real Imag 

K = 

= 0.2 



-.0961 

-.4811 

-.0528 

-.4773 

-.0837 

-.5270 

-.6108 

-.3625 

-.5422 

-.3721 

-.6270 

-.3965 

-.1247 

-.1151 

-.1176 

-.1210 

-.1270 

-.1266 


K = 0.3 


-.0690 

-.7736 

-.0489 

-.7260 

-.6471 

-.5592 

-.5015 

-.5691 

-.1344 

-.1729 

-.1168 

-.1820 


-.1211 


















4.3 GRADIENT OF STEADY MEAN MASS FLUX VECTOR 

Before proceeding to a discussion of results from thick wing or wing-body flow problems, a 
discussion of one of the influences of the underlying steady mean component of flow is 
required. An additional approximation is introduced by deleting one of the boundary condition 
terms which is not computed with sufficient accuracy by the panel method of ref. 1 . 

The aerodynamic surface and wake surface boundary conditions shown by eq. (31) contain 
the following term involving the gradient of the mass flux vector: 

(d • V w,) • [ D -V (v0s --1 . 

The value of this term depends on the second order derivatives of steady mean velocity 
potential, and its accurate evaluation places stringent requirements on the accuracy of the 
steady mean flow solution. As noted in the preceding, the steady mean component of flow 
is computed by the panel method of ref. 1 . If that method is not sufficiently accurate, 
then this term cannot be included in the unsteady panel method formulation. 

The accuracy of the steady panel method was tested by computing second order derivatives 
of the velocity potential at the surface of a sphere in an uniform, incompressible flow. 

These computed values were then compared with values computed from the exact solution 
to this flow problem, i.e.. 



where a is the radius of the sphere and r is the distance from the center of the sphere to 
the point where the potential is evaluated, viz., 

r = (x2 + y2 + z2)‘'^. 

The spherical surface was approximated by 968 panels. Twenty-two panels were uniformly spaced 
in rows along meridian lines running from the forward stagnation point to the aft stagnation point 
and there were forty-four such panel rows around the circumference of the sphere. The problem 
was solved using Morino type boundary condition, cf. sec. 3.4; hence, the potential and its 
first and second order derivatives with respect to the local panel coordinates were readily 
evaluated at each panel center from the solution doublet distribution. 

The exact solution to the flow problem was differentiated with respect to the compressibility 
coordinates, and the components of the derivatives were transformed to the local panel 
coordinates. The derivatives of the potential, with respect to the local panel surface coordinates, 
should then compare with those computed by the panel method. For a typical meridian 
line row of panels the results are as shown in table 7. 


63 


Table 7. — Computed Vs. Exact Values of Velocity Potential 
and Its Derivatives at Surface of A Sphere 


0(deg.l* 

0 


0 

n 

hi 


0 

W 

4.09 

.49835** 

.02183 

-.00004 

-.84453 

-.03974 

.06799 

.50130*** 

.03542 

.00007 

-1.51171 

.00001 

-1.51163 

12.27 

.48831 

.10618 

-.00271 

-.47719 

.01194 

-.87561 

.49126 

.10565 

.00018 

-1.48240 

.00004 

-1.48164 

20.45 

.46835 

.17415 

-.00406 

-.46848 

.01045 

-.81333 

.47134 

.17400 

.00026 

-1.42409 

.00010 

-1.42214 

28.67 

.43886 

.23910 

-.00530 

-.43757 

,00909 

-.76325 

.44191 

.23920 

.00031 

-1.33751 

.00017 

-1.33408 

36.82 

.40041 

.29923 

-.00601 

-.38881 

.00749 

-.69249 

.40350 

.29997 

.00033 

-1.22384 

.00024 

-1.21893 

45.00 

.35381 

.35313 

-.00632 

-.33226 

.00580 

-.61314 

.35683 

.35505 

.00032 

-1.08480 

.00031 

-1.07870 

53.18 

.29998 

.30278 

.40006 

,40317 

-.00606 

.00029 

-.26680 

-.92269 

.00413 

.00038 

-.51912 

-.91597 

61.36 

.23999 

.24242 

.43905 

.44327 

-.00519 

.00024 

-.19756 

-.74049 

.00460 

.00045 

-.41078 

-.73386 

69.55 

.17510 

.46891 

-.00415 

-.13673 

,00543 

-.30231 

.17697 

.47433 

.00018 

-.54180 

.00049 

-.53603 

77.73 

.10662 

.48936 

-.00291 

-.08068 

.00555 

-.19175 

.10778 

.49553 

.0001 1 

-.33082 

.00053 

-.32659 

85.91 

.03591 

.49993 

-.00129 

-.02526 

.00543 

-.07203 

.03631 

.50633 

.00004 

-.11223 

.00055 

-.11003 


* Meridian angle from forward stagnation point 
** Computed value 
*** Exact value 


Except at the first panel, which is closest to the stagnation point, the value of the potential and 
its gradient are computed with sufficient accuracy to be regarded as a satisfactory steady flow 
solution for the purpose of evaluating the surface pressure distribution. The second order 
derivatives, however, are considerably in error over the entire meridian line. Because of these 
large errors, the boundary condition term described at the beginning of this subsection was 
deleted in the boundary conditions of the demonstration cases. 
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4.4 APPLICATION OF THE PANEL METHOD TO 
WINGS OF FINITE THICKNESS 

Figure 34 shows the panel arrangement of a test case wing having finite thickness. This 
wing has the planform shown by fig. 18 and an airfoil section that is 4.8 percent thick at 
a point 40 percent of the chord aft of the leading edge. The thickness of this wing is 
sufficiently small that the lifting pressure induced on its surface as a result of harmonic 
flow incidence is accurately predicted by the thin wing theory method of ref. 8. Using 
the paneling shown by fig. 34 the lifting pressure has been computed, satisfying the unsteady 
thick wing flow theory shown by eqs. (30) through (35). The wing is oscillating in pitch 
about a spanwise line through its mid-root chord point with a reduced frequency of 0.3577. 

The lifting pressure coefficient distributions along four chord lines are shown by figs. 35 
through 38. These figures provide a comparison of the real and imaginary parts of the 
pressure with those computed by the method of ref. 10. Since the thick wing panel method 
should yield results in this case which are nearly identical with results from the method of 
ref. 8, the close correlations shown by fig. 35 through 38 imply validation of the panel 
method. 



Figure 34. — Pane! Arrangement for Thick Test Wing 
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Figure 35. - Thick Wing Unsteady Lifting Pressure Comparison at 18. 1% Semi-span 
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Figure 36. — Thick Wing Unsteady Lifting Pressure Comparison at 5 1.2% Semi-span 
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Figure 37. — Thick Wing Unsteady Lifting Pressure Comparison at 81.7% Semi-span 
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Figure 38, — Thick Wing Unsteady Lifting Pressure Comparison at 97,7% Semi-span 



This wing flow problem was computed using the Morino-type boundary conditions described in 
one of the examples of sec. 3.5. The wing is represented by distributions of unsteady sources 
and doublets on the upper and lower surfaces of the wing and on a surface which closes its tip. The 
strengths of the sources are set equal to the unsteady flow incidence and the unsteady disturbance 
velocity potential is required to vanish at all points interior to the wing surface. A direct applica- 
tion of these boundary conditions in the panel method, however, does not lead to an accurate 
representation of the flow in the neighborhood of the wing trailing edge. This problem was 
corrected by a modification of the Morino-type boundary conditions; the nature of the problem 
and its remedy are described in the following. 

The panel method is applied by defining three networks of panels covering the wing surface: 

(1) the upper surface, (2) the lower surface, and (3) the tip closure surface. As indicated 
above these are combination source-doublet networks; and, in the terminology of ref. 2, 
sec. 2.4, these three networks are typed as follows: 

NTS(K) = 1 and NTD(K) = 12 for K= 1,2,3. 

These networks have control points at the panel centers and at the mid-points of the panel edges 
which are also network edges, cf., fig. 15. Two boundary conditions are required at each panel 
center control point and a single boundary condition is required at each network edge control 
point. 

Again, using the terminology of ref. 2, sec. 2.4., the Morino-type boundary conditions are imposed 
by setting the boundary condition parameters as follows: 

at panel center control points (ref. 2, tables 2.2 and 2.3): 

NLOPTl = 5 
NROPTl = 5 
NLOPT2 = 7 
NROPT2 = 2 

at network edge control points: 

NLOPTl = 0 
NLOPT2 = 7 
NROPT2 = 2 

This specification of the boundary condition parameters causes the source strength to be set equal 
to the flow incidence at panel centers and causes the disturbance potential to vanish at limiting 
points approaching the control points from the interior of the wing. 

The doublet strength values at control points at the edges of two adjoining networks are used, 
in part, to enforce continuity of the doublet strength between the adjoining networks across their 
abutting edges, ref. 1, sec. 6.0. In the present case of the thick wing, the doublet strength on the 
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wake along its upstream edge is forced to match the sum of the doublet strengths of the upper 
and lower wing surface networks along their line of abutment with the wake. At the edge control 
points of the upper and^lower surface networks, however, if the Morino-type boundary conditions 
are imposed, causing to vanish at these points, the Kutta condition fails to be correctly satis- 
fied. This failure is a failure of the boundary conditions to correctly determine the direction of the 
disturbed flow at the wing trailing edge. The flow should be tangent to the wake surface at this 
edge; but, because eq. (82) is not satisfied precisely at the wake line of abutment, the flow direction 
can be in a direction significantly different from that of the tangent to the wake surface. This 
failure is remedied by requiring the flow to be tangent to the upper surface of the wing at the edge 
control points. 

At the edge control points along the trailing edge of the upper surface network, the boundary con- 
dition parameters (ref. 2, tables 2.2 and 2.3) are given the following values: 

NLOPTl =0 
NLOPT2=3 
NROPT2 = 2 


As a result, 

A 

wg • n = 0 

at edge control points of the upper surface along the wake and this condition causes the unsteady 
flow to have a direction tangent to this surface along its trailing edge. 
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4.5 APPLICATION OF THE PANEL METHOD TO A 
WING-BODY-TAIL^NACELLE CONFIGURATION 


Figure 39 shows the paneling for a typical subsonic, twin engine transport with the following 
characteristics: 

Wing area = 297.0 square meters (2759 square feet) 
aspect ratio = 8.71 
taper ratio = 0.267 
sweep @ C/4 = 31.5 degrees 
dihedral = 6 degrees 
root chord = 856.7 cm (337.3 inches) 
tip chord = 228.6 cm (90.0 inches) 

MAC = 603.18 cm (237.474 inches) 



Figure 39. — Paneling of Typical Transport Configuration 
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Horizontal tail area = 89.3 square meters (830 square feet) 
aspect ratio = 4.50 
taper ratio = .250 
sweep @ C/4 = 33 degrees 
dihedral = 7 degrees 
root chord = 662.4 cm (260.8 inches) 
tip chord = 165.6 cm (65.2 inches) 

MAC = 463.6 cm (182.5 inches) 
span = 1 863.0 cm (733.4 inches) 
tail arm = 2053 cm (808.4 inches) 
tail volume coefficient = 0.926 

Body cross section maximum width = 502.9 cm (198.0 inches) 
maximum height =541.0 cm (213.0 inches) 

The body is represented by sources and doublets distributed directly on the actual surface 
except on the closure section aft ot the horizontal tail. A cylindrical wake surface with gen- 
erator parallel with the undisturbed freestream emanates from the aft body surface simulating 
flow separation from the body surface at the body station where the trailing edge of the 
horizontal tail intersects the body. The body closure surface is simply a flat surface normal 
to the treestream and at the leading edge of the body separation wake. Morino-type boundary 
conditions, cf. sec. 3.4, are applied at the body surface; hence, the source distribution strength 
is set equal to the complex amplitude of the unsteady flow incidence at the body surface and 
the unsteady potential is set to zero at the body interior. Source and doublet panels are also 
placed on a surface which closes the aft body interior to the cylindrical wake separation surface. 
This body closure surface is a flat surface normal to the freestream direction. The unsteady 
disturbance potential is set to zero at either side of the body closure surface, thereby, causing 
the doublet strength to vanish but providing a source distribution together with control points 
along the edge of the body closure surface adjacent to the wake surface. These boundary 

conditions lead to undisturbed freestream flow interior to the cylindrical wake separation 
surface. 

At the wing, tail, and nacelle surfaces thin lifting surface boundary conditions are applied. 

The effects ot thickness and the effect of a mean steady component of flow are ignored. 

A nacelle is treated, essentially, as a ring wing, and it is on a strut having a wake. 

Figure 40 shows the paneling of the transport used in a doublet-lattice computation. The 
body paneling shows the surface of interference for the body. The actual body was represented 
by unsteady line doublets to generate its slender body theory aerodynamic influence. 

The transport is assumed to be oscillating in pitch about a point 848.3 inches aft of the nose 
of the body. The reduced frequency was chosen as 0.0 1 74 so as to be typical of the short 
period frequency ot an aircraft of this type. The flow Mach number was chosen as 0.87. 

The complex amplitude of the lift and pitching moment coefficients induced by the pitch 
oscillation were computed by both the panel method and the double-lattice method. The 

results are summarized as follows: 

panel method doublet-lattice method 

(3.5033) + i(. 0647) 1/2 = (3.0022) + i(. 1301) 

l/2Cjyj= -(15.4742) -i(.75279) l/2cJi= -(13.2033) -i(.9273) 

Where the pitching moment is about a point 234.95 cm (92.5 inches) forward of the nose. 

The solution to this problem required approximately 1 000 CPU seconds on the Ames CDC 
7600 computer. 






5.0 COMPUTER CODE USER’S GUIDE 


The computer program of the unsteady panel method is a modified form of the steady flow panel 
method computer program described by ref. 1 and 2.f Input data for otecuting the unsteady flow 
program is nearly identical with the input data required by the steady flow program; hence, 
this section is written as a supplement to the user’s manual of the steady flow program, viz., 
ref. 2. The reader is expected to have ref. 2 in his possession and to be familiar with its contents. 

There are three principal differences between the input requirements of the two computer pro- 
grams. The differences consist of the following additional data required by the unsteady flow 
program : 

(1) boundary condition data describing unsteady surface motion at the control points 

(2) wake paneling grid point data providing for harmonic variation in the transverse 
component of wake vorticity, and 

(3) mean steady component of flow data at the control points of the unsteady flow 
problem. 

The input data is described in detail in sec. 5.2. 

Output data is also similar to the output data of the steady flow program. The primary difference 
is that the output data of the unsteady flow program is complex, i.e., the program generates 
the real and imaginary parts of the complex amplitudes of the unsteady flow parameters. 

The output data is described in sec. 5.3. 

5.1 GLOBAL COORDINATE SYSTEM 

All geometric input describing network grid points is expressed in the global coordinate system, 
see sec. 1.0, ref. 2. The global and compressibility coordinate systems are related by a user 
supplied angle of attack (ALPC) in degrees and angle of sideslip (BETC) in degrees. From 
these two angular rotations, the program computes the components of a unit vector in the 
direction of the x-axis of the compressibility coordinate system. This unit vector is as follows; 

c = COMPD(l)t+ COMPD(2)| + COMPD(3)k (88) 

AAA 

where i, j, k are the unit base vectors of the global coordinate system. 

5.2 INPUT DATA REQUIREMENTS 
5.2.1 STEADY MEAN COMPONENT OF FLOW DATA 

By reference to the boundary conditions equations, i.e., eqs. (3 1 ) and (34), it is seen that 
at each control point the following steady mean component of flow data is required: 


For computer system requirements refer to sec. 1.5 of ref. 2. 
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steady mean density, 


9 


M. 


r- 


u. 


A . 

C • 


(89) 


and steady mean mass flux vector, 

= • (90) 

In addition, the gradient of the mass flux vector is required to evaluate the second term 
of the boundary condition shown by eq. (31 ). That term is evaluated from second order deriva- 
tives of the steady velocity potential, with respect the local panel coordinates, see app. E. 

The user may elect to ignore the influence of the mean steady component of flow. This 
choice is set by the logical variable LSTDY in subroutine INPUT. If the statement 
LSTDY=.TRUE. appears in subroutine INPUT, the program expects mean steady flow data. 

If the statement LSTDY=.FALSE. appears, the steady mean density and steady mean 

mass flux vector are given their freestream values. The default value of LSTDY is .FALSE. 

The steady mean flow data is provided by executing the steady flow program of refs. 1 and 2. 

The steady flow program must be executed using precisely the same geometry and paneling as 
that to be used in the unsteady flow program. The control point numbering and the control 
point locations are therefore identical with those to be used in the unsteady flow problem, and 
the following data is saved for each control pointf: 

JC, control point number; 

ICP, panel number containing JCth control point; 

PLCL(I), 1=1, 2, 3, local panel coordinates of the control point; 

TSC( 1 ), source strength ; 

TSC(2), doublet strength; 

TSC(I), 1=3, 4,5 components of the gradient of the doublet strength in 
the global coordinate system; 

W(l), average potential; 

W(I), 1=2, 3,4 components of the perturbation mass flux vector in the global 
coordinate system; 

DDO(I), 1=1,2 components of the gradient of doublet strength in local 
panel coordinate system; 

X1(I),I=1,2,3 second order derivatives of the doublet strength with respect to the 
local panel coordinate system. 

In the unsteady flow program these quantities for each control point are contained in * 

labeled common block STDY. 


t This data must be saved in the subroutine OUTPUT, see fig. 1 . 1 2 of ref. 2, using the 
following statements: 

For each case: 

WRITE (8,8001) AMACH, ALPHA (lACASE), BETA (lACASE) 

For each panel center control point: 

WRITE (8,8001) JC,ICP.PLCL,TSC,W,DDO, (XI (I), 1=4,6) 8001 FORMAT (4020) 
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The data is contained on an octal formatted file identified to the prograni as tape 13. The first 
record on the file contains the following values: AMACH,ALPHA(IACASE), BETA(IACASE). 
Subsequent records contain the data which is listed above and which appears in labeled 
common block, STDY. Tape 13 is read as a sequential file in the unsteady flow program and it 
is generated by the steady flow panel method program by inserting the previously mentioned WRITE 
statements. The steady flow program must be executed with NEXDGN=1 and bit number one must 
be set in the variable IFLAG (e.g., IFLAG = IFLAG. or .1). 

5.2.2 BOUNDARY CONDITIONS 

Boundary conditions are specified at each control point following the procedures of sec 2.4 
of ref. 2. The general form of the boundary condition is identical with that of ref. 2, except 
that the flow parameters are complex as shown by eq. (79), i.e., 


Cu(wu* • ns) + f^^ + + cj(wj2* •rig)+fg-v0g* + Djj0g* = BET.* (91) 


In certain cases the coefficient BET* is represented by the value of the mass flux boundary 
condition shown by eq. (31). In these cases the coefficient BET* is computed by Subroutine 
CBET by choosing option 5 of table 2.3, ref. 2. As shown by eq. (31), the normal component 
of mass flux is equated to 


BET*=[-W^.(0*X^^^)-(d*vW 3) -n^ + PsUj^*] (92) 


where R is the position vector of the control point in the global coordinate system. This 
equation contains the followingjcomplex amplitudes of the^ unsteady surface motion 
parametei|: surface rotation, *; surface displacement, D*; and surface rate of displace- 
ment, Uj^ . These parameters are related to time dependent quantities as follows: 

R (e* , D = R (d* ei^O ^ and u„ = R (u„* e‘^t) 

Where co is the circular frequency in radians per second and R( ) denotes the real part. 
The circular frequency is input to the program in the following ratio: 


OMGBAR = w/Uq = d3 


(93) 
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which is related to the reduced frequency as 

k = c3c/2 (94) 

where c is a characteristic length expressed in the units of length used in defining the 
freestream velocity magnitude. From eqs. (93) and (94) it is clear that the input parameter 
is related to the reduced frequency of the unsteady motion as follows: 

OMGBAR = 2 k/c . (95) 

As an example of the computations carried out by subroutine CBET, consider the case 
of the configuration shown by fig. 41 undergoing pitch oscillations about its center of 
mass, i.e., its eg. 



Figure 41. — Airplane Configuration Oscillating in Pitch 
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Assume the following values: 

Xq = 9.14 meters (30 ft), 

Uq = 213.4 m/sec (700 ft/sec), 

1^0 1 = 3° 

C = 7.92 meters (26 ft), and 

k - 1.0 

and assume that the x axis is parallel with the freestream and that the x axis passes through 
the eg. The frequency input value is given by eq. (95) as follows: 

OMGBAR =2.0(1.0)/13.0 
= 0.07692 


The surface displacement at the control point with position R is 



D*(R)- 0 oX(R-Ro) 

(96) 

where 

^0 = (0.0)t+ (3.0/57.296)t+ (0.0)1^ 
= (0.0)t+ (0.05236)? + (0.0)k 

(97) 

and 

Ro=(30.0)t+(0.0)?+(0.0)i^. 

(98) 

The surface rotation at the control point is 



e (R) = 0 o 

(99) 

The surface displacement rate at the control point is 



u„* = id5UoD*(R)-n3. 

(100) 

Substituting into eq. 92 and using eq. 89 and 90, the right hand member of the boundary 
condition is seen to be given by 

r* 


BET* = 

- ^ + V0g - $ Mq2 $ . V 03 ) . X ft J 


- 




/ \ 1 

^ *^<^s) Uq (^ 0 * ^ (R - Ro)) • "s e-‘"C‘R 

(101) 
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This expression is evaluated in subroutine CBET at each control point storing the complex 
amplitude BET* in the labeled common block BCOND. 

Evaluation of BET* for the example requires the FORTRAN code shown by fig. 42; however, 
the second term contained in the brackets of eq. (101) is not included in this computation. 

As shown by sec. 4.4, second order derivatives of the velocity potential are not evaluated with 
sufficient accuracy for this term to be included in the computations, and it is not included 
in any of the examples computed for this report. 

As seen in fig. 42, the factor 


g-iw^RMQ2/|32 (102) 

is computed by a call to subroutine CMAB, computing the scalar product 

(FACTOR)=c-R 

where the components of c are contained in COMPD(I) while the components of R 
are contained in ZC(I). The system subroutine CEXPl then computes the factor shown by 
eq (102). The subroutine CBET next computes the following cross products: 

D = 9oX(R-Ro) 

^xv ^ 

and the subroutine CRZAB computes the scalar products: 

CFAC = c • X rig) and DN = D • rig . 

Finally, the operation 

CSUM = FSVM(I)*(DN-DFAC), 

appearing just before the text concerning the logical variable LSTDY, is equivalent to 

CSUM = X (R - Rq)) • rig - c • (^t*X tig)j . (103) 

This expression represents the bracketed quantity of eq. (101) when the influence of the steady 
component of flow is ignored. 
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SUBROUTINE CBET 

COMIVION/ACASE/ALPHA(5),BETA(5),FSVIVI(5),FSV(3,5),IACASE,NACASE 
COMPLEX THETA 
COMMON/BCD AT/R (3) ,THETA(3, 5) 
COMMON/CNTRO/ZC(3),ZNC(3),ZDC,IPC,ICC,JZC,JCN,KC 

common/comprs/amach,betams,betam,sbetam,abetms,alpc,betc, 

COMPD(3),AROTC(9),AROTCI(9),CZINV(3,3) 

COMPLEX BET 

C0MM0N/BC0ND/CU,CL,TU(3),TL(3),DU,DL,BET(5),NCT,NL0PT,NR0PT,NBIN 

common/freqdt/omgbar,omegb,omg 

LOGICAL LSTDY 

COMMON/OUTDAT/IFLAG,LSTDY 

COMMON/STDY/JC;ICP,PLCL(3),TSC(6),W(4),DDO(2),X1(3) 

COMPLEX CFAC,CSUM,CXV,CEXPI,CZEXP,D,DN 
DIMENSION CXV(3),D(3),PHES(3),IMP{3) 

DATA THETA/(0.,0.), (.05236, 0.), (0., 0.)/ 

DATA R/30.,0.,0./ 

CALL CMAB(ZC,C0MPD, FACTOR, 1,3,1) 
CZEXP=CEXPI(FACT0R*0MEGB*(BETAMS-1.)) 

GA=CU+CL 

GD=.5*{CU-CL) 

DO 1000 l=1,NACASE 

D(1)=THETA(2,l)*(ZC(3)-R(3))-THETA(3,l)*(ZC(2)-R(2)) 

D(2)=THETA(3,I)*(ZC(1)-R(1))-THETA(1,I)*(ZC(3)-R(3)) 

D(3)=THETA(1,I)*(ZC(2)-R(2))-THETA(2,I)*(ZC(1)-R{1)) 

CXV(1)=THETA(2,I)*ZNC(3)-THETA(3,I)*ZNC(2) 

CXV(2)=THETA(3,I)*ZNC(1)-THETA(1,I)*ZNC(3) 

CXV(3)=THETA(1,I)*ZNC(2)-THETA(2,1)*ZNC(1) 

CALL CRZAB(D,ZNC,DN,1,3,1) 

DN=DN*CMPLX(0.,OMGBAR) 

CALL CRZAB(CXV,C0MPD,CFAC,1,3,1) 

CSUM=FSVM(I)*(DN-CFAC) 

IF(.N0T.LSTDY)G0 TO 1000 

CALL CSCAL(BETAMS,C0MPD,ZNC,TMP) 

CALL CMAB(ZNC,TMP,FACT0R, 1,3,1) 

CALL VADD(TSC(3),TSC(1 )/FACTOR,ZNC,TMP,3) 

DO100J=1,3 

100 PHES(J)=GA*W(J+1)+GD*TMP(J) 

CALL CMAB{PHES,COMPD,FACTOR, 1,3,1) 

FACTOR=FACTOR*(SETAMS-1.) 

CALL VADD(PHES,FACTOR,COMPD,TMP,3) 

CALL CRZAB{CXV,TMP,CFAC, 1,3,1) 

CSUM=CSUM+DN*FACTOR-CFAC 
1000 BET(I)=CSUM*CZEXP 

RETURN 
END 


Figure 42. — Example Subroutine CBET 
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The logical variable LSTDY is supplied to the program as one of the input parameters in 
subroutine INPUT and it is contained in the common block labeled OUTDAT. If this variable 
has the value TRUE, then subroutine CBET computes the steady mass flux vector and the 
steady disturbance density ; it then forms the following result: 


CSUM= - (Uq c+V(^3-cMo2g X rig) 

* ® i“ U„( 9„*X (r - Ro))- fis 

This operation is expressed as 

CSUM = CSUM -(v0g - c Mq2 $ • v0g) • (d* X n^) 

-Mo^c • v0giw(^oX (R-Rq)) • tig 


(104) 


( 105 ) 


where CSUM on the right is given by eq. ( 103). Referring to figure 42, the call to subroutine 
CSCAL leads to an evaluation of the conormal to the surface; hence, the components of TMP(3) 
contain the components of n^ after the call to CSCAL. The call to CMAB leads to an evaluation 
of the scalar product of the surface normal with the surface conormal, i.e.. 


(FACTOR) = n • n^ . (106) 

The call to VADD leads to an evaluation of v 0g at the upper side of the surface. Recalling 
eq. (74), viz., 

'^u==(^)a+^("c‘' + ^^) (107) 

and noting that 


TSC (2 + J) = VMs for J = 1 , 2 , 3 

while 


TSC(l) = ag , 


the call to VADD is seen to yield 


"c ''s +^^s ; 


further, since 


CA = 1.0 and CD = 0.5, 

regardless of whether Morino-type boundary conditions (CU=1.0, CL=-1.0) or mass flux boundary 
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condition (CU-1.0, CL=0.0), the “do loop” leads to 

PHES(J)=v0g 

evaluated at the upper side of the surface. The next subroutine call to CMAB 
evaluates the following scalar product: 

C * V0g 

and the next operation yields 


FACTOR = FACTOR*(BETAMS-1.0) 


In turn, the calls to VADD and CRZAB yield 

(108) 


TMP (I) = v0g - c c • 

(109) 

Recalling that 

CFAe=(v^,-eM„2«.vO-(9o’x«s)' 

(110) 

the final operation, viz., 

DN = S,.(9*X(R-R„))i5, 

(111) 


CSUM - CSUM+DN* FACTOR - CFAC, 

(112) 


yields an evaluation of eq. ( 1 05 ). 

The final operation of subroutine CBET, yields the boundary condition value, viz., 

BET* = CSUM(e~i^c-RMo2/^2y 


The preceding example describes a single boundary condition consisting of the rigid body 
mode shape: pitch oscillation. Provision is made for up to five boundary condition cases 
representing five different mode shapes; thus, for example, there is provision for introducing 
five different centers of rotation for the pitch oscillation or for introducing five completely 
different mode shapes. One must recognize, however, that the flow Mach number and reduced 
frequency are fixed and that these parameters are the same for each boundary condition case. 

Two labeled common blocks, ACASE and BCON, contain the information related to multiple 
cases. ACASE contains the variable NAGASE defining the total number of cases to be executed; 
it also contains the variable lACASE identifying the case currently being executed and three 
scalars and a vector which can be evaluated for each case. BCON contains the boundary condi- 
tion defining parameters associated with any control point which are the coefficients appearing 
in eq. (91) and up to five values for BET* are allowed for each of the two boundary conditions 
at a control point. 
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5.2.3 UNSTEADY WAKE REPRESENTATION 

The unsteady wake representation available within the program provides capability for construe? 
ting an elaborate model of the unsteady wake vorticity convection. Wake surfaces are represented 
by doublet networks which, as noted in sec. 3.2, can be made to form a stream surface of the 
steady mean component of flow. The unsteady vorticity can be convected along the mean steady 
streamlines on the wake surface with the average mass flux of the mean steady flow. 

As shown by eq. (37), viz.. 


T(iw/^2'J J d|/wj 

jLi* (C , s) = ju* (c , s) e c 

the wake model determines the doublet distribution on an entire wake network in terms 
of the doublet strength along the upstream edge of the network. The integral contained 
in eq. (37) yields the time period, i.e., 

1 £ 

At= - f d^/w, (114) 

t 

representing the period of time required for unsteady vorticity to be convected from point c 
on the network edge to the point £ interior to the network. Writing eq. (37) in terms of the 
time delay, viz., 

ju* (C , s) = M* (c , s) [cos (ojAt) - i sin (wAt)] , 

the doublet distribution is seen to be a harmonic function of the time period. Since the 
steady disturbance velocity is at least an order of magnitude less than the freestream velocity, 
the factor Wg is nearly unity (cf. the discussion preceding eq. (37) in sec. 3.1 ); hence, the 
unsteady wake vorticity varies nearly as a harmonic function of the x compressibility 
system coordinate, i.e., 

jLt* (x , s) ^ ju* (c , s) j^cos (x - c)) - i sin (co (x - C))] (115) 

where 


X-C = Uq At . 


The doublet strength distribution on the wake is represented in terms of doublet panel 
distributions. As noted in sec. 3.4, a doublet panel distribution is a quadratic function 
of the local panel coordinates, c.f. eq. (50), and the quadratic is determined by its least 
squares fit to the values of the actual doublet distribution at the evaluation points shown 
in fig. 1 2. As shown above, at a wake surface the actual doublet distribution is very nearly 
a harmonic function of the downstream distance, see fig. 43. A quadratic can approximate a 
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Lifting surface Wake surface 



Figure 43. — Wake Panel Spacing 


harmonic function to within three significant figures over a domain which is one sixteenth 
of the wave length of the harmonic function ; therefore, the length of a wake panel in the 
streamwise direction should not be greater than one sixteenth of ( 27 r/oo) in regions where 
that level of accuracy is required. 

The wake networks are defined in subroutine INPUT and they are paneled like any other doublet 
network defining the grid points of the panel corner points, sec. 1.1, ref. 2. Control points are 
defined and boundary conditions are specified only along the upstream edge of a wake network, 
sec. 1.3.2 of ref. 2. The time period contained in eq. (37), forevaluation the downstream wake 
doublet strength, is introduced at each doublet singularity parameter point in subroutine DASPL. 
Eq. (37) is expressed as 

M* (1) = M* (J)e'K<^/^^)TDLY(I) 

where 

C(l) 

TDLY(I)= f d^/w^ (117) 

m 
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where C(J) is the wake streamline coordinate of the network edge singularity parameter upstream 
of the Ith singularity parameter on the wake, the latter, having the wake streamline coordinate 
fi(I), fig44. 



Figure 44. — Singularity Parameters Along Wake Streamline Coordinates 

In its present form, subroutine DASPL generates TDLY(l) as follows; 

TDLY (1) = c • (r( 1) - to) (118) 

where c is the unit vector in the compressibility direction, eq. (88), while R(I) and R(J) are 
the position vectors, respectively, of the downstream and upstream singularity parameters. 

This expression approximates the wake vorticity convection time period, ignoring the effect 
of the average mean steady disturbance mass flux vector. This approximation is used because 
the analysis capability of the steady flow panel method does not currently provide for an 
evaluation of eq. (108). 

5.3 OUTPUT DATA DESCRIPTION 

With the exception of pressure data, the output data of the unsteady flow panel method is a subset 
of the steady flow panel method as listed in Table 3~1 of ref. 2. Of course, values of all unsteady flow 
parameters are expressed as the real and imaginary parts of their complex amplitudes. The following 
pressure coefficients are output. 
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Linear, Unsteady Disturbance Pressure 



Upper Surface 

Lower Surface 

Difference 

Real 

CPLINU 

CPLINL 

CPLIND 

Imaginary 

CPLNUI 

CPLNLI 

CPLNDI 


C * 





Quadratic, Unsteady Disturbance Pressure 



Upper Surface 

Lower Surface 

Difference 

Real 

Imaginary 

CP2NDU 

CP2DUI 

CP2NDL 

CP2DLI 

CP2NDD 

CP2DDI 


2 

= (iwps0VPo+(Ws/PoUo)‘V'^*) 

where 





Linear, Total Disturbance Pressure 



Upper Surface 

Lower Surface 

Difference 

Real 

Imaginary 

CPLNUT 

CPLUTI 

CPLNLT 

CPLLTI 

CPLNDT 

CPLDTI 


Cp*=Cp=«= (119)- -^ 6' 

where Cp* (119) denotes the value of Cp* computed from eq. (119) 

Quadratic, Total Disturbance Pressure 



Upper Surface 

Lower Surface 

Difference 

Real 

Imaginary 

CP2DUT 

CP2UTI 

CP2DLT 

CP2LTI 

CP2DDT 

CP2DTI 


(119) 

( 120 ) 

(17) 

( 10 ) 

( 121 ) 
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Cp* = Cp* ( 1 20) - ^^^2 J + iJ »s 


(122) 


whe re 


^ 9 A ^ . 


( 11 ) 


and Cp* (120) denotes the value of Cp* computed from eq. (120) 
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6.0 CONCLUSIONS 


An advanced panel method for solving subsonic unsteady potential flow boundary value prob- 
lems has been presented. The method is shown to be applicable to small disturbance flows when 
the flow can be described with adequate accuracy by a solution to the linear, first order 
approximation to the unsteady potential flow equation. The method is also shown to be 
applicable to arbitrary configurations. 

The principal feature of the method which required validation is a kernel function approximation. 
Solutions to the unsteady flow equation are related to solutions to Helmholtz’s equation 
by a coordinate transformation and a change of dependent variable. Using this relationship, 
a solution to an unsteady flow problem is expressed in terms of an integral equation derived 
from Helmholtz’s theorem. The panel method is based on a reduction of this integral 
equation to a system of linear, complex algebraic equations and the reduction is accomplished 
by a typical panel method approximation leading to a panel by panel evaluation of the integral 
expression appearing in the integral equation. The kernel function approximation is required 
for the closed form evaluation of the integrals. 

The validity of the kernel function approximation has been examined both theoretically and 
empirically. The theoretical examination consists of a direct error analysis. The empirical 
examination consists of comparisons of computed panel method results with results from 
alternative computational methods. The alternative methods do not contain kernel function 
approximations of the same nature as that of the panel method. The comparisons include 
solutions to unsteady flow problems of thin, planar wings and for T-tails; they demonstrate 
that the panel method solution tends to converge to the correct solution as the number of 
panels is increased. 

In the analysis of unsteady flow on thin, planar wings the panel method yields aerodynamic 
forces which converge nearly monotonically with increasing panel density. Monotonic con- 
vergence for a sequence of thin wing solutions tended to occur only for a particular set of 
panel density distributions. In the case of the more complex T-tail, however, the panel method 
tended to converge with increased panel density but convergence was not always monotonic. 

This condition is thought to result from the particular choice of panel density distribution. 
Finding the panel density distribution leading to monotonic convergence may not be plausible 
except for very simple configurations undergoing very simple modes of motion. 

The validity of the panel method for computing unsteady flow about bodies of finite 
thickness was tested by an indirect comparison. The unsteady flow was evaluated about a 
wing having a finite (but small) thickness (viz., 4.8 percent) with the wing oscillating in pitch. 

The unsteady flow was evaluated for the case of thick wing theory with the panels located on 
the actual wing surface and with the flow boundary conditions satisfied at the actual wing 
surface. The unsteady flow was evaluated a second time but to provide a comparison case 
using thin wing theory. In thin wing theory the panels lie on a flat mean wing surface 


and the flow boundary conditions at the actual surface are satisfied only to a first order 
approximation; but, since accurate alternate methods are available for its evaluation, thin 
wing theory provides a reliable comparison case. The lifting pressure distributions for 
these two cases are compared and found to be in very close agreement. This close agreement 
shows that the panel method is valid for bodies of finite thickness. 

Application of the method to complex configurations was demonstrated by computing the 
unsteady flow about a wing-body-tail-nacelle configuration undergoing pitch oscillation. 
Results from this computation are compared with equivalent results computed using the 
well known doublet-lattice method. The comparisons are in sufficient agreement to indicate 
validity of the panel method and the capability to analyze very complicated geometries. 
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APPENDIX A 

GENERAL THEORY OF FLOW 


A.l INTRODUCTION 

This appendix contains a derivation of the general theory of compressible, inviscid, nonconducting, 
irrotational, isentropic flow expressed as a boundary value problem. The derivation leads to a 
single, nonlinear partial differential equation governing the flow; and it leads to the analytical 
expressions for boundary conditions appropriate to solid, aerodynamic surfaces and to wake 
surfaces. The derivation proceeds from results appearing in ref. 14 and its purpose is to 
provide a convenient reference for the development of the approximate theory of flow appearing 
in app. B. 

A.2 FLOW EQUATION 

The derived governing equation (called, as usual, the flow equation) is a single equation in 
terms of a single flow parameter, viz., the velocity potential, but it is derived here from a 
system of six equations in terms of six variables; the components of velocity, pressure, density, 
and internal energy. The system of six equations, taken as the primitive statement of the 
theory, represents the conservation laws (in their forms appropriate to an inviscid, nonconducting 
fluid undergoing an isentropic process) along with the constitutive relation for the fluid, the 
latter being an equation relating internal energy to the pressure (a mechanical variable) and the 
density (a deformation variable). The flow equation is derived by a process of reduction wherein 
the primitive system of equations are combined so as to eliminate the explicit appearance 
of the original variables. In the reduction process two additional variables are introduced by 
two additional equations which express the constraints imposed to obtain irrotational and 
isentropic flow. 

In the statement of the primitive system of equations and in the derivation involving those 
equations, the laws of conservation of mass, momentum, and energy are stated independently 
of one another, for example, the law of conservation of momentum is expressed without 
asserting that mass is conserved. This independence is retained because of the particular 
requirements of the development in app. B which leads to an approximate theory. In the 
development contained in app. B approximations are postulated regarding the conservation 
laws. The interdependence of those approximations is examined using the equations of app. A; 
therefore, the laws of conservation are initially stated (cf., eqs. (A. 2a) and (A.3a)) as inde- 
pendent relations. 


The primitive system of equations is as follows: 


• Continuity equation (eq. 5.4), ref. 14— 


This equation follows from the law of conservation of mass and relates the dilitation 
rate of the fluid to its rate of change of mass density as 


where 


v-v 

f> dt 


d 9 ^ ^ 

— = — + V„ • V. 
dt 9t ^ 


(A. la) 


Euler’s equations of motion eq. (6.9), ref. 14- 

“Vc 1^ n 

— + — vp = 0. 

dt p 

These equations follow from the law of conservation of linear momentum when 
mass is conserved causing the law of conservation of linear momentum, viz., 


(A.2) 




+ p — +vp = 0, 
at 


(A.2a) 


to reduce to eq. (A.2). 

Energy equation eqs. (34.2) and (34.3), ref. 14- 


P— +PV* V =0. 
dt ^ 


(A.3) 


This equation follows from the law of conservation of energy when eqs. (A.l) and (A.2) 
are satisfied causing the law of conservation of energy, viz. , 


® ^ ^c • ^c) • ^c) + ^ ■ 


c 1 _ \ de 


— + — vp 1 + p — + p V V = 0 , (A.3a) 
dt p / dt 


to reduce to eq. (A.3). 


• Constitutive relation 


de_ 1 d 
dt 7-1 dt 


(P/P) 


This equation is derived from the perfect gas relations (sec. 31, ref. 14) viz. 

dt '' dt 


(A.5) 
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where 


R = Cp - Cy and 7 = Gp / Cy . 

Eqs.(A. 1 ), (A.2), (A.3), and (A.4) are a determinate system of six equations in terms of 
six flow parameters, namely, 

P,P,e,V^. 

The specific internal energy e is eliminated as an explicit dependent variable by combining 
eqs. (A.3) and (A.4) to obtain 

Id 

— ; — (p/p) = -(p/p)v- V . 

7 -I dt 

Introducing eq. (A. la) to eliminate the dilitation rate, after some manipulation, the pressure 
and density are found to be related as follows: 

dp _ dp 

dt dt 

This result is integrated, following the motion of a specific but arbitrarily chosen fluid 
particle, to find the well known isentropic pressure-density relation, viz., 


(P/^o) = (P/Po)^- (A. 8 ) 

Equations (A.l), (A.2) and (A.8) form a determinate system of five equations in terms of five 
flow parameters, namely, 

P , P , Vj, . 

For the case of irrotational flow the system of equations governing the flow may be combined 
to eliminate explicit appearance of the velocity by introducing a new flow parameter, namely, 
the velocity potential. Letting the flow be a disturbance from a umform freestream and 
introducing the kinematical constraint of irrotationality, viz., vX = 0, the flow velocity 
may be expressed in terms of a potential as follows: 


Also, as shown by eqs. (3,5) and (17.1), ref. 14, the fluid particle acceleration may be expressed as 

dV^ 

=v\h + Zo X V„ 

dt 

where 

-j._-j.-k 90 I -1. 

co = vx Vg and 0 = — +-V^ • ; 

thus, tor irrotational flow in which o3 vanishes, the acceleration terms of Euler’s equations of 
motion, viz,, eq. (A. 2), can be expressed as the gradient of a scalar: 

J^c • Vc) + -vp = 0. (A.IO) 

Noting, further, that eq. (A.8) supplies a unique relation between pressure and density 
(i.e., the flow is barotropic), the following relation follows from the rule for differentiating 
a definite integral taken along an arbitrary path in space from a point where the pressure 
and density are uniform: 


/ 


p(x,y,z,t) 


dp 

P(P) 


p(x,y,z,t) 

/ ^ 


1' 


P(P) 


dp + 


1 


p(p(x,y,z,t)) 


vp(x,y,z,t) 


1 -k 1 -k 

^Po=T^P 


P(Po) ^ ^ 


wherevp(p) vanishes because the dependence of pressure on density is independent of 
position and v p^ vanishes because the lower limit of the integral is assumed to be at 
a spatial point where the pressure is uniform. Using this result, Euler’s equations of 
motion are expressed entirely as the gradient of a scalar quantity: 







Integrating this expression over an arbitrary path from a point in the undisturbed fluid, 
Bernoulli’s integral to Euler’s equations of motion is found by introducing eq. (A.9) and 
choosing the potential 0 to be zero in the undisturbed fluid. The result appears as 
follows: 


D0 

Dt 




V0 + 


/dp 
/ — = 0 
Po " 


(A.ll) 
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where 


= 

Dt at o 


’ V0 . 


Further, substituting eq. (A.9) into the continuity equation yields the following expression 
in which the velocity components are eliminated as explicit variables: 



dt ■ 


Finally, eqs. (A.8), (A.l 1), and (A.l 2) provide a determinate system of three equations 
in terms of three flow variables, namely, 

P , P , 0 • 


The pressure and density are eliminated as explicit variables by introducing the speed of 
sound as a new variable (c.f. sec. 35, ref. 14) viz., 

2 dp 
c = — ■ . 
dp 

Combining eq. (A.7) with eq. (A. 1 3), the speed of sound is found to be given by the 
following ratio: 

c^ = 7 p/p . 


Since eq. (A. 14) is a property of the fluid, the differential with respect to time following 
a fluid particle can be expressed as 


dt p dt 


and, on integrating, following the motion of a specific but arbitrarily chosen particle, 
the speed of sound is evaluated as follows: 



c^ - Cq2 = (7-1) 



P 


where and appear as a result of having chosen a path of integration beginning at 
an undisturbed point in the flow. This result is combined with eq. (A.l 1) to obtain the 
speed of sound in terms of the velocity potential, viz., 


c2 = Cq2 


/D0 1 _ 

(T-1)( — +-V0 • v0 


(A. 12) 


(A. 13) 


(A. 14) 


(A.15) 
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Pressure and density are eliminated by using eq. (A. 13) to express the continuity equation 
(as given by eq. (A. 1 2)) in terms of the speed of sound. The result is as follows: 



1 dp 
p dt 


Differentiating Bernoulli’s integral (as given by eq. (A.l 1)) and combining the result with 
this form of the continuity equation, a determinate system of two equations, consisting 
ofeq. (A.l 5) and 


1_ d/D0 

2 dtW 


+ — V0 • V0 


')■ 


is found in terms of two flow parameters, namely, 


<j) and c. 


(A. 16) 


The desired result, viz., a single equation in terms of a single flow parameter, is obtained 
by eliminating the speed of sound between eq. (A.l 5) and eq. (A. 16) to find the following 
equation in terms of the disturbance velocity potential : 


9 I 7 / D0 1 . . 

- (7-1) {— +- • V0 


-1 d /D0 1 _ _ 


‘) 


this equation is termed the “flow equation.” 


(A. 17) 


The fluid pressure is evaluated in terms of the velocity potential by substituting the 
isentropic pressure-density relation (i.e., eq. (A.7)) into Bernoulli’s integral (i.e., eq. (A. 1 1 )) 
and carrying out the indicated integration. When the result is expressed in terms of a 
pressure coefficient, it appears as follows; 


Cp=(p-Po)/ 



2 




7/7-1 


- 1 . 
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A.3 BOUNDARY CONDITIONS 

The equations of the preceding section govern the kinematics and dynamics of the fluid 
flow; as such, they represent conditions which must be satisfied at every point of the flow 
field. This section introduces additional conditions which must be satisfied on surfaces 
which form the boundaries of the flow field. 

The bounding surfaces are of three types. They consist of the aerodynamic surface of the 
aircraft, the surfaces of its wake, and a far-field surface which completely surrounds the 
aerodynamic surface, but which is at a large distance from it. The aerodynamic surface has 
the characteristics of the surface of a solid body; therefore, fluid particles do not penetrate 
this surface. In addition, the aerodynamic force on an aircraft is the result of the fluid 
pressure acting on the aerodynamic surface. The wake has two surfaces which are adjacent 
to one another and are separated by an infinitesimal distance. These wake surfaces form 
two sides of a vortex sheet; fluid particles do not pass through them and the component 
of velocity tangent to them is discontinuous. Disturbances to the fluid flow are required 
to propagate toward the far-field surface and are required to become vanishingly small 
there. 


A.4 AERODYNAMIC SURFACE 

Let S denote the aerodynamic surface of the aircraft at the present instant of time, and 
denote the volume enclosed by this surface as V. In the analytical representation of the 
flow problem the fluid is continuous and fills the entire space; hence, the surface S 
encloses the mass of fluid, M, computed as follows: 


M = 


fffp dv 


V 


In the case of an arbitrary motion of S, the mass M may vary with time; and, assuming 
mass is conserved, this variation implies that there is motion of fluid through S in violation 
of the aerodynamic surface boundary condition. Spatial points of the surface S have a velocity 
denoted as Tl, see sec. 3. 1 ; hence, the time rate of change of the mass contained in V is given by 


^ _ £im 
dt "At 


^ (x,y,z,t + At) dv + JJp (x,y,z,t +At) At Ij • h ds v,z,t)dvj| 

v' S' v' 

Jffjt 


s' 


where S' and V' are spatially fixed but coincident with the moving surface and volume (S and V) 
at the instant of time t, see fig. 45. 
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Figure 45. — Aerodynamic Surface Variation in Time 

Since mass is conserved, the continuity equation, i.e., eq. (A.l), is introduced and the volume 
integral contained in eq. (A. 19) is as follows: 

XCT " -fSf^ ' (f> ^c)dv . 

V V 

Applying the divergence theorem leads to 

///I7 dv--)y*pV^ -nds 
V' S' 


so that on substituting this expression into eq. (A.l 9), the time rate of change of mass is 
found as follows: 

• (A.20) 

This equation shows that the time rate of change of enclosed mass is evaluated by the 
surface integral sum of 



If there is to be no exchange of fluid from V, interior to S, with the fluid exterior to S, 
then this differential quantity must vanish at every point of S; therefore, the aerodynamic 
surface boundary condition is given by 

pu^ = M • n on S (A.21) 

where 



A 

u * n 
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and 




M=pV^ (A.22) 

is the fluid mass-flux vector. 

In the case of steady flow, wherein u^j = 0, the boundary condition is found by 
integrating the steady flow continuity equation over V to find 

///v-Mdv = 0. 

V 

Applying the divergence theorem, the aerodynamic surface boundary condition is found 
to be 

M • n = 0 on S . (A.23) 

The aerodynamic force acting on the surface S, see fig. 46, is found as 

Fs = -//pnds. (A.24) 

S 

Let S(, be a control surface, moving with the fluid and surrounding S, and let V be the volume 
of fluid between and S. The total force acting on this volume of fluid must vanish ; hence, 

F = -iy*pnds- // pn ds + ^ /ffo Vgdv =0 (A.25) 

^c S ^ V 


Applying the transport theorem (ref. 14, sec. 4), the volume integral becomes 

Iff yy*pV^Vc*nds+ //pVgV^.*nds 


dt 


V 


Sc 


where S[. is chosen to be a spatially fixed surface coinciding with at time t. In steady flow 
the first two terms on the right vanish (the first by definition and the second by virtue of 
eq. (A.23)); hence, 

— ///p dv = //p • n ds 

V Sc 


Substituting this result and eq. (A.24) into eq. (A.25), the force acting on S is found in 
terms of an integral over the spatially fixed control surface, , as follows; 


Fs = 


= // (p Vc • n + pn 


S' 


^ds . 


(A.26) 
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V 



Figure 46. — Control Surface for Computing Aerodynamic Force 

Hq. (A. 26) also follows directly from Euler’s equations of motion and the continuity 
equation. In steady flow these expressions are given by 

Mv-V^,+ vp = 0 (A.27) 

and 

V • M - 0 ; (A.28) 

thus, eq. (A.27) can be written as follows: 

V • (m Vj.) +vp = 0 . (A.29) 

Integrating over the volume V and applying the divergence theorem, eq. (A.29) leads to 

77* |^(Vc m) • n + pnjds = 0 . 

S,+S 

Introducing the steady aerodynamic surface condition shown by eq. (A. 23) and the definition 
shown by eq. (A. 24), the aerodynamic force on S is found again as 

Ps= //(^ » + 

In steady flow, therefore, if the mass flux vector, velocity vector, and pressure satisfy eqs. 

(A.27) and (A.28), then eqs. (A. 23) and (A. 26) are valid expressions for the aerodynamic 
surface conditions. 


The above results are well known and appear in most standard texts on fluid mechanics. 
They are developed here to provide a special emphasis for app. B, where approximations 
are introduced. The approximations are imposed on the relationships governing the pressure 
and density in the flow. In steady flow they lead to an approximate mass-flux vector which 
satisfies eqs. (A.28) and (A.29) exactly even though mass and momentum are not conserved 
exactly. As a result, the surface conditions shown by eqs. (A. 23) and (A. 26) are taken to 
be valid in the approximate theory of flow. 
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A.5 WAKE SURFACE 


The wake surface is a vortex sheet containing the vorticity shed from the trailing edge of a 
lifting surface; and as such, it is the idealization of the real flow which was originally 
introduced by Prandtl (ref. 1 5 sec. 1 9. 1 ). The boundary layers at either side of a lifting 
surface separate from the surface and form a thin wake region of viscous flow trailing 
downstream. Like the thickness of the boundary layers, the thickness of the thin wake 
region of viscous flow tends to zero in the limit of zero viscosity. As noted by sec. 2.7, 
ref. 7,. there is also a jump in potential across the wake surface along its upstream edge 
where it is attached to the lifting surface, and this jump in potential is equal to the circu- 
lation on the lifting surface. Since, in general, there can be both a spanwise and timewise 
variation in circulation on the lifting surface, there can be a gradient in the jump in potential 
along the wake surface. This gradient represents the vorticity shed from the lifting surface 
into the wake surface. 

The wake is a free vortex sheet except along its upstream edge where it is attached to the 
lifting surface; also, the wake surface does not propagate through the fluid with at least 
the speed of sound. It follows from these two conditions (ref. 14, sec. F) that the fluid 

pressure is continuous across the wake surface, i.e., 

1 p 1 = 0 on W. (A.30) 

Additionally, in view of the discontinuity in potential across the wake, there may be a 
discontinuity in fluid velocity across the wake, viz., 

I ]]:^0 onW. (A.31) 

Eqs. (A.30) and (A.31) cause the wake to be a particular type of surface of discontinuity. 

Formal procedures for analyzing the flow in the neighborhood of a surface of discontinuity 
are described by ref. 16, chapter C, and those procedures make use of the geometry shown 
by fig. 47. In fig. 47 the quantity S denotes a material surface surrounding the material 
volume V. This volume of fluid is divided into two parts, V*" and V", by^the surface 
s and this dividing surface need not be a material surface. The volume V is surrounded by the 
surface S’*'+s and the volume V" is surrounded by the surface S" + s. Finally, the speed of 
surface displacement along the surface normal is given by the following expression : 



Figure 47. — Surface of Discontinuity 


u+-ft = 


A 

u • n 


V^ - n on S+ 


Uj^ on s 


• n on S" 


Uj^ on s . 


As shown by ref. 14, eq. (6.1), momentum is conserved interior to the surface S if 


— Sffp Vc dv - - JJp ft ds. 


Introducing the transport theorem (ref. 14, eq. (4.1)), the volume integral can be expressed 
as follows: 


^jfffoVc<lv = 


///[^ 

\ 7 I- 


s(pVc) 


V V 

and, by applying the divergence theorem, it follows that 


+ V^ • v(pvJ+pV^ V. 


dv; 


dt' 


ff, 

V V s 

Introducing the volume subdivision, these expressions yield 

ff. 


P • n ds . 


V 


dt 


V 


: at 




V,V,vftds 


+ .^OV^V^-ftds-//(pV^)\„ds+ fffvfv,. 


ds 


ff P ^ ■ ff P n ds ; 

C+ (- 1 - 
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hence, letting and S“ approach s and noting that and V“ vanish in the limit, 
while noting the arbitrariness in the choice of the integral over s, this equation becomes 

^c(^c * " “ %)1 + IIPI n = 0 . (A.32) 


A development identical to that above, but applied to conservation of mass, leads to the 
Stokes-Christoffel condition (ref. 16, eq. 189.14), viz., 

|'p(Vc-fi-uJ|-0. (A.33) 

This result may be used in conjunction with eq. (A.32) to obtain 

IPJ n-p*U± [vj = 0 (A.34) 

where U is the speed of propagation of s: 


Introducing the requirement that the pressure be continuous across the wake, 
viz., eq. (A. 30), it follows from eq. (A.34) that 

P±U±[V^.J=0 on W. 

From eqs. (A.34) and (A. 30) it follows that only the tangential component of velocity 
can be discontinuous across the wake surface, viz., 

[v^ • n]-0 onW 

and 

[ • t 0 on W 

where t is an arbitrary unit vector tangent to the wake surface. These equations also show 
that the discontinuity in tangential component of flow must vanish unless the speed of 
propagation vanishes at either side of the wake surface ; hence, 

p— Uj^ = • n on W, (A. 35) 

Eq. (A. 35) is an analytical statement of the requirement which must be satisfied for the wake 
surface to be a material surface. 
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Let Bernoum’s integral, i.e., eq. (A.l 1), be evaluated at either side of the wake surface and 
form the difference of the two resulting integrals to obtain 

p" EpI 

• |v0] + [ v0 • v0] = - J dp/p + y* dp/p = - y* dp/p 

o o o 


where and p” denote the pressure at either side of the wake surface. As a consequence 
of eq. (A.30), it follows that 



and, recognizing that 





• V0 + — V0 • V0 


Bernoulli’s integral is found to lead to 
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v^ll=0 on W 


(A36) 


This condition must be satisfied if the wake surface is to be, simultaneously, a material 
surface and a surface of velocity discontinuity. 


In the following the wake boundary condition shown by eq. (A. 36) is developed into a partial 
differential equation governing the value of [[01 on the wake surface. This development 
follows a consideration of the kinematics of the fluid flow at the wake. To this end con- 
sider the two positions of the wake surface shown by fig. 48. One position is that of the wake 
surface at time t; the other is its position a short interval of time, At, later. The two 
positions are denoted W(t) and W(t + At). 


q' 



Figure 48. — Wake Surface Motion 
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The unit vector, n, normal to W(t) at Q is directed at a surface point Q' of W (t + At). 
Let Q have the compressibility system coordinates (x,y,z) and let the compressibility 
system coordinates of Q' be denoted as (x + Ax, y + Ay, z + Az). The quantities 
(Ax,Ay,Az) represent the coordinate differences of the two surface points. 

Let A|[0] denote the difference in the value of 10]] at points Q and Q' so that 

AI01 ^ 11 01' - m 

At At 

At 

^ 0'+-(^+ 0'~-0~ 

At " At 

-[A0/At] 

where the prime denotes evaluation at the point Q' while the absence of the prime denotes 
evaluation at the surface point Q. Taking the limit as At tends to zero leads to the result 

where the operator 6 / 5t denotes the time rate of change apparent to an observer moving 
with the component of velocity of the surface along its normal. 

With 0 expressed as a function of the compressibility coordinates, the change in the 
value of 0 during the interval of time At is given by 

A 90 

A(^) = At Uj^ n • v0 + ^ At ; 

hence, dividing by At and taking the limit as At tends to zero results in 

60 ^ ^ 90 

- = u„n-v0 + — • 


Evaluating eq. (A. 38) at either side of the wake surface and forming the difference of the 
two expressions yields 



A 


n 



Introducing eq. (A. 37) into eq. (A. 39) leads to 



(A.37) 


(A.38) 


(A.39) 


(A.40) 
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This result is called the kinematical condition of compatibility (ref. 16) and -this condition 
must be satisfied if the discontinuity is to persist in time. 

Noting that 

n *[v0]i =|n • vj 

the requirement that the normal component of velocity be continuous 
leads to 

n • [v^l = 0 


and the kinematical condition of compatibility given by eq. (A.40) reduces to 



on the wake surface. The wake boundary condition given by eq. (A.36), therefore, can 
be expressed as 

4 +i[v^ . V J = 0 on W , 

The objective, here, is to reduce the wake boundary condition to a partial differential 
equation in terms of the discontinuity in the velocity potential. The kinematical condition 
of compatibility shown by eq. (A.42) has accomplished this_ objective for the first term. 
The second term is expressed in the desired form by expressing |Vc J in terms of the 
average flow velocity at the wake surface, viz., 




Using this definition, deviation velocities at the wake may be defined as 

q+ = V/-V,^ and q-- 

with the result that 

■q''' = -q‘ 

The deviation velocities contain the wake velocity discontinuity; hence, the average flow 
velocity is continuous at the wake, i.e., 

[vJ = 0onW; 


(A.41) 


(A.42) 


(A.43) 


(A.44) 


(A.45) 


(A.46) 


(A.47) 


106 


and it is tangent to the wake surface, i.e., 

• n = 0 on W. 


Consider the second term in the wake boundary condition as shown by eq. (A.43) 
apply the above conditions. Since 

M = V -Vc'- V cn W. 

introducing eq. (A.44) leads to 

M = . 2V„ • . 5- . ^ 

and, because eq. (A.45) leads to 



it follows that 

M = 2V^-Iqlon W. 


Further, since 

Uq+V 0- Vj^ + q , 

then 

|q]=|v0|. 

It, therefore, follows that 

= 2 M ; 

and, on applying Hadamard’s lemma given by pp. 492 and 493, ref. 16, this result 
reduces to 


(A.48) 


(A.49) 


Finally, combining this result with eq. (A.43), the wake surface boundary condition is 
found as 

7 : 1*^1 ‘ = 0 on W . 

ot 

This is the desired result: a partial differential equation governing l<j)l onW. 


The wake surface is shown by eq. (A.5 1) to be a surface of convection on which the 
quantity |[0]1 varies in time as a result of [[01 being convected downstream with the 
mean flow velocity tangent to the wake surface, W. Given the wake location, the 
average flow velocity, and the value of 101 along the upstream edge of the wake, eq. 

(A.5 1 ) can be integrated to find the value of [[01 at any point on the surface. The upstream 
edge of the wake surface is the trailing edge of a lifting surface. The spanwise value of 
101 along this edge is known to be equal to the spanwise value of the circulation on the 
lifting surface (ref. 7, eq. 7-39). Eq. (A. 51), therefore, does not represent an inde- 
pendent boundary condition. The flow conditions on the wake are completely determined 
by eq. (A.51) in conjunction with the flow equation and the independent boundary condi- 
tions on the aerodynamic surface. 


APPENDIX B 

APPROXIMATE LINEAR THEORY OF FLOW 
B.l INTRODUCTION 


An approximate, linear theory of flow is formulated in this appendix starting from the funda- 
mental postulate that all nonlinear terms of the flow equation, viz., eq.(A. 17), are sufficiently 
small as to be negligible. This postulate reduces the flow equation to the following linear 
approximation: 


9 1 D^<p 

20 = 


Co 2 Dt 


2 ' 


(B.l) 


Recall that the flow equation (for the case of irrotational flow) is equivalent to the 
following system of six equations (cf, app. A): 


continuity, 




(A. la) 


conservation of linear momentum. 


dVc 


+ p — +vp = 0; 


conservation of energy, 




de 

+ p — + pv* V_ = 0 ; 
dt c 


and the ideal gas relation, 


de 1 d 

T," 7. ^P/P^ ■ 

dt 7-1 dt 


(A. 2a) 


(A.3a) 


(A.7) 
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Since the fundamental postulate is imposed on the exact flow equation, the approximations 

to each of the equations in this equivalent system of six equations are indeterminate; 

also, since both the aerodynamic surface and wake surface boundary conditions follow 

from the requirements of the laws of conservation of mass and linear momentum at those surfaces 

(cf. sec. A.4 and A.5), the approximate flow theory cannot be derived without choosing 

additional postulates which determine the approximations to those two conservation laws.f 

The additional postulates are chosen assuming the unsteady flow to have a steady mean 
component so that the flow may be viewed as an unsteady small disturbance superimposed 
on a steady small disturbance flow. The additional postulates are chosen such that the 
steady mean component of the flow satisfies the boundary value problem of the linear, 
steady flow theory of ref. 1, summarized in sec. B.2. As shown in sec. B.2, the additional 
postulates are stated as approximations to two relations: (1) the relationship between mass 
density and velocity potential, viz., 

p = pW , 

and (2) the relationship between pressure and velocity potential, viz., 

P = P(0) • 

These approximate relations are chosen such that mass and linear momentum are conserved 
when there is no unsteady component of flow. 

Consistent with these objectives, the unsteady flow is described by a velocity potential such 
that the unsteady flow velocity is given by 

= +% + ♦„) (B.2) 

where is the velocity of the uniform freestream, 0^ is the velocity potential of the 
unsteady small disturbance component of flow, and </>g is the velocity potential of the 
steady small disturbance component of flow. 

B.2 MEAN, STEADY COMPONENT OF FLOW 

For the steady mean component of flow to be predicted by the theory of ref. 1 it must satisfy 
the following boundary value problem derived in sec. 3 of ref. 1 : 

Governing differential equation — 

V • Wg = 0 

where 

W^^Po(Uo + W3) 


(B.3) 

(B.4) 


+ This is a significant point and should be recognized as such. The additional postulates 
lead to the so-called “mass flux boundary conditions” and the use of the second order 
approximation to Bernoulli’s integral to Euler’s equations of motion. 
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and 



As shown by ref. 1 , when (j>^ satisfies eqs. (B.5), (B.6), (B.8), and (B.9), the steady mean flow 
has the following characteristics: 

( 1 ) Mass is conserved because the continuity equation, (A. 1 a), is satisfied. 

(2) Conservation of linear momentum, eq. (A. 2a), reduces to Euler’s equations of 
motion, viz., 

^Ps " ® > (B.IO) 

because the continuity equation is satisfied. 


Ill 


(3) Euler’s equations of motion are satisfied when the pressure coefficient is 
approximated by eq. (B.8);and, therefore, Euler’s momentum theorem holds, cf. 

ref.14. 

(4) A solution to the linear flow equation (viz. eq. (B.5)) in the_domain V surrounded 
by the closed surface E is uniquely determined if 0g and Wg are uniquely 
specified on E. 

Several additional characteristics of the flow are important to the purpose of the following 
section, i.e., the derivation of the unsteady flow theory. One of these characteristics is the 
approximation to the relation between mass density and velocity potential. That approxi- 
mation is derived by expressing the continuity equation as 


v^. VP3. (B.ll) 

^S 

This expression shows (cf. ref. 14, sec. 3 and 4) that the steady dilitation rate is proportional 
to the rate of convection of the fluid mass density. The steady dilitation rate is also evaluated 
by writing the flow equation (i.e., eq. (B.4)) as follows: 


These two expressions for dilitation rate combine to express the rate of mass density con- 
vection as follows: 


7 VPs ih) 

^s 


XX 


(B.12) 


This result can be expressed in terms of the disturbance mass density, viz., 

^s ^s “ ’ 


(B.13) 


and that expression appears as 

=r('-"s7Po+---)(Uo+^*s)- ’"s' 

where the series is the result of a binomial expansion. The order of magnitude of the 
disturbance mass density is evaluated from this expression to find 

Ps'/"o'=''(<'sAo)' 

thus, to a first order approximation 
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Writing eq. (B.14) as 



integrating with respect to x, and evaluating the constant of integration in the freestream, 
the desired approximation to the relation between the flow mass density and velocity is 
found to be as follows: 

Mq2 


Substituting this approximation into the expression for mass flux, viz., 




and neglecting terms of second order in leads to the approximate mass flux vector 
shown by eq. (B.4); hence, eq. (B. 1 5) is consistent with the theory of ref. 1 . 


The derivation of the unsteady flow theory also requires an evaluation of the order of 
magnitude of the disturbance pressure. That evaluation follows from the law of con- 
servation of energy (viz., eq. (A. 3a)) and the ideal gas relation (viz., eq.(A.4)). Since both 
mass and linear momentum are conserved, conservation of energy (i.e., eq. (A.3a)) reduces 
to the usual energy equation, viz.. 


As a consequence, the energy equation can be combined with the ideal gas relation to 
find the isentropic pressure-density relation, cf. app. A., 


(Ps/Po) = (Ps/PoF • 


The isentropic pressure-density relation is now expanded in a Taylor series about the 
freestream conditions to find 


Ps"Pq 

Po 




+ ^(7-1) Co 


2 


Ps-Pq \^ 

^o / 


+ ••• 


where 

Cq2 =(dp/dp)o = 7(Po/Po)- 


(B.15) 


(B.16) 


(A.8) 
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On neglecting terms which are shown by eq. (B.15) to be second and higher order in 
a first order approximation's found, viz,, 


Ps-Pq 

Pq 



/Ps “P o 



(B.17) 


Combining eq. (B.l 7) with eq. (B.l 5), the order of magnitude of the disturbance pressure 
is found to be given by 


Ps “ Pq Pq ‘"o ^ (^s/^o) • 


(B.l 8) 


B.3 UNSTEADY FLOW APPROXIMATIONS 


Approximations for the unsteady flow theory, which are consistent with the approxima- 
tions of the preceding, are deduced by following the procedure leading toeqs. (B.15) and 
(B. 18). The continuity equation is expressed in terms of the divergence of the mass flux, i.e.. 



(A.l) 


and an unsteady disturbance mass flux is deduced by writing the mass flux in terms of 
steady and unsteady components, viz.. 


p Vc = Ws + w^, , 


(B.l 9) 


so that the continuity equation becomes 


ap 


(B.20) 


because the steady flow component satisfies the continuity equation, i.e., 

W^^O. 


The unsteady flow equation, i.e., eq. (B.l), is now written as follows: 




Uq Dt 


it(DtW)- 


Comparing this form with eq.(B.20), the following hypothesis is asserted: 


M 


o D 


Wu-Poiv^u - i-0-^ 57 (M 


(B.21) 
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On integrating eq. (B.22) with respect to time, the unsteady mass density is found to be 
given by the following: 


P 




+ P 


s ■ 


As a test of the validity of the above hypothesis, consider <j>^ to be independent of 
time. The disturbance mass flux approximation shown by eq. (B.21) then reduces to that 
of the steady flow theory, viz., that shown by eq. (B.4);also, the disturbance mass density 
approximation shown by eq. (B.23) reduces to that of the steady flow theory shown by 
eq. (B.15). The hypothesis therefore is seen to be consistent with the steady flow theory 
of ref. 1 . 


For the purpose of deducing an approximation to tlie unsteady pressure coefficient, 
recall Bernoulli’s integral to Filler’s equations of motion, viz., 

P 

D0 1 /'dp 

+-v0*v</)+ / — = 0. 

Dt 2 J p 

Po 

Assuming that Filler s equations of motion are satisfied at least to a first order approxi- 
mation, the isentropic pressure-density relation can be expressed as follows: 




P -Pf 


+ - (7-1) Cq 


2 /P Po 


Substituting eq. (B.24) into eq. (A. 1 1 ), the integral with respect to pressure is evaluated 
to find 



(B.22) 


(B.23) 


(A. 11) 


(B.24) 


(B.25) 
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Bernoulli’s integral, in terms of this series, is inverted by letting 


P-Pq 

Po 


oo 


n=l 


a„ (W))" 


where 


D0 1 ^ 

f(0) =— 


and evaluating the coefficients using the uniqueness property of the power series. 
The result of substituting those coefficients into eq. (B.26) provides the desired result, 
namely. 


P-Pq 

Po 




1 ^ _ 
— + — ' v<j> 

2c^-\Dt 2 




Introducing 0 = ^u neglecting terms which are of second order in , the 
unsteady disturbance pressure is found by subtracting eq. (B.8). The result is the 
following approximation to Bernoulli’s integral; 

Mo^ \ ^ \ 

which, by reference to eqs. (B.4) and (B. 1 5), can be expressed as 

p - Ps « - ^Ps (-^u)^ + • ^(<^u)) • 

Comparing eq. (B.28) with the usual approximation of linear, unsteady flow theory, viz., 

P - Ps ^ - (Po + Po Uq (^u)^)’ 

the freestream values of the mass density and the mass flux are seen to be replaced in the 
present theory by the values of those quantities evaluated in the steady mean component 
of flow. As a point of information which should be recognized, eq. (B.28) represents an 
approximation to Bernoulli’s integral to Euler’s equations of motion. In the present 
context Euler’s equations of motion are approximated. Equation (B.28) therefore is the 
appropriate approximation for evaluating aerodynamic forces at the boundaries of a flow; 
however, if the pressure is viewed as one of the flow parameters, then the pressure is 
related to the velocity potential by eq. (B. 24) truncated as 


P-Ps 

Po 


(B.26) 


(B.27) 


(B.28) 


(B.29) 



and eq. (B.23) expressed in terms of the total disturbance density, i.e., 

p-pQ D0 

————— ii. . 

Pq ■ Uq 2 Dt 

Combining these two expressions leads to the expression for the unsteady pressure 
evaluated as a parameter of the flow, viz., 


(B.30) 


P - Po D0 

This result is identical with eq. (B.29) (i.e., the usual approximation of linear, unsteady 
flow theory), and this is the expected result in view of the fundamental postulate of the 
approximate theory, namely, that terms, having an order of magnitude smaller than 
O(0/Uo^’ neglected in the governing differential equations. Because of the funda- 
mental postulate, Euler’s equations of motion are represented only to a linear, first 
order approximation in (0 / Ug); products of 0^ and 0^^ appear in eq. (B. 28) only 
because eq. (B.28) is an integral to a linear system of equations in which 0g appears as a 
parameter. 


B.4 AERODYNAMIC SURFACE BOUNDARY CONDITION APPROXIMATION 

The aerodynamic surface boundary condition, which was derived in app. A, states the 
following requirement: Fluid particles located at the aerodynamic surface S must have a 
motion satisfying eq. (A. 21 ), viz. 


pUj^ = W • n . 


(A.21) 


An approximate form of this equation is now derived for the case of unsteady, small 
disturbances from a steady mean flow. 


The unsteady flow is a consequence of a small, unsteady disturbance motion of the 
surface S^relative to its steady mean location , fig. 49. A point P on S has 
position R and becomes the point Q on with position when the unsteady 
motion ceases. The point Q is said to be undergoing a displacement, i.e., a change 
in position, computed as follows: 


D = R-Rs. 


(B.31) 
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The displacement of points on Sg are assumed to be sufficiently small that the value of a 
flow parameter evaluated at P can be approximated by a truncated Taylor series expansioi 
about Q, e.g., , 


(j) (P, t) « 0 (Q, t) +(yd>)Q • D . 



Figure 49. — Unsteady Surface Displacement 


The velocity of the point P relative to the point Q is given by 

^_3D 

“~at 

and the unsteadiness of the motion is assumed to be such that its Fourier representa- 
tion has a reduced frequency band whose upper bound has the order of magnitude 
unity or less. As a consequence, the velocity of a surface point, when normalized by 
the freestream velocity, has an order of magnitude which is equal to or less than that 
of its displacement. Further, because of the assumed small displacement, the rotation 
of the surface at point P may be linearly approximated, viz., 

^«-vXD; 

2 

and the unit vector normal to S may be approximated, viz.. 




(B.32) 


(B.33) 


(B.34) 


(B.35) 
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In turn, the normal component of surface velocity is approximated as follows: 

A 3D ^ 0D 
= n • ■— • — . 


Substituting eqs. (B.35) and (B.36) into eq. (A.21 ), the surface boundary condition is 
approximated as follows: 

p Uj^ « W • + 0 X ng^ on S . 

Introducing eqs. (B.l 9) and (B.23), this result becomes 

{'>u+Ps)i'„^(Ws + 5u)-(Ss + «Xn^)on S 


% 


where 




^ D<t> 
2 ^ 


u 


The truncated expansion of the flow parameters, of, eq. (B.32), leads to an expression 
evaluated on the steady mean surface, viz., 

^s^n^Wg + e X -vWs) • ng + Wy • n^ on 

where all nonlinear terms involving unsteady quantities are neglected as small; further, as a 
consequence of eq. (B.6), the unsteady disturbance mass flux is found to be required to 
satisfy the following boundary condition on the steady mean surface: 

Wu • ng^-Ws -(0 X ns)- (D • vwj • ng + pgu„ on 


where 




^s (Uq ^ (^s) X 



(B.36) 


(B.37) 


(B.38) 

(B.21) 

(B.4) 

(B.15) 
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B.5 WAKE SURFACE BOUNDARY CONDITION APPROXIMATION 

The jump in potential across the wake surface is represented by a doublet sheet; hence, the 
wake boundary condition, i.e., eq. (A. 51), is expressed as follows: 

8ju 

- — + *\7jji-0 on W 

St ^ 

where 

is the strength of the doublet sheet. If eq. (B.39) is multiplied by the flow mass density, 
the boundary condition is an expression in terms of mass flux, i.e., 

Sfx ^ 

P — +W,^ -VM = 0 on W . 

5t ^ 

This result follows from continuity of pressure across the wake surface, i.e. 

II p I = 0 on W, 

and the isentropic pressure-density relation, i.e., eq. (A.8); whereby, 

I p n = 0 on W 

so that the average mass density across W is equal to the mass density at either side 
of the wake. 


The flow parameters appearing in eq. (B.41) and evaluated at a point P on W can be 
expressed by a Taylor series expansion about the steady mean location of P on W^. 
The result appears as follows: 

/ - ^ \/3m 3D ^ - -.3/Lt \ 

(p + D-vp+--)I- + --v« + D-v- + ---j 

^ D W|^ + '••)• V + D ■ v/i +•••) = 0 on . 

Introducing . 

P = Pu + Ps > = (W + ( Wy)^ , and w = p, + p„ 


(B.39) 


(B.40) 


(B.41) 


(A.30) 
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and noting that jUg satisfies the steady mean wake boundary condition, viz., eq. (B.7), after 
deleting products of unsteady quantities, the unsteady wake boundary condition is 
expressed as follows: 


dt 



= 0 on Wg. 


(B.42) 


Recalling eq. (B.40) and comparing eq. (B.42) with eq. (B.28), the unsteady wake boundary 
condition, as approximated by eq. (B.42), imposes the requirement that the approximate 
disturbance pressure be continuous across the wake. The approximation shown by 
eq. (B.42) is therefore consistent with the approximation for pressure. 
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APPENDIX C 

DERIVATION OF INTEGRAL EQUATION FORMULATION OF 
THE UNSTEADY FLOW PROBLEM 


The objectives of this appendix are twofold. The first objective is to provide a derivation 
of the integral equation which is shown as eq. (38). The second objective is to express that 
integral equation in terms of panel integrals by introducing a change of variable and a 
coordinate transformation. These panel integrals are evaluated by the methods of app. D 
to generate the unsteady aerodynamic influence coefficients which appear in sec. 3.4 

In sec. C. 1 the unsteady flow equation, shown as eq. (30), is shown to be identical in form 
with a particular form of Helmholtz’s equation. This form of Helmholtz’s equation follows 
from the application of the Prandtl-Glauert transformation to Helmholtz’s equation when it is 
expressed in rectangular Cartesian coordinates. The appUcation is such that two units of 
measure for distance are introduced into the formulation. One of these is the measure 
associated with the rectangular Cartesian coordinates of Helmholtz’s equation; the second 
is that of the compressibility coordinate system in which the flow problem is defined, t 

Section C.2 introduces the theorems of Helmholtz which provide an integral equation 
satisfying Helmholtz’s equation. The integral equation is expressed in invariant form 
such that it can be expanded in terms of any coordinate system related to a rectangular 
Cartesian reference frame by an admissible transformation of coordinates. The Prandtl- 
Glauert transformation is admissible; and, in sec. C.3, the form of the integral equation 
satisfying the flow equation is derived. This derivation includes the definition of 
the conormal of the surface of integration. 

Section C.4 delineates the line vortex integrals which arise at the perimeter of a surface 
having a doublet distribution. 

Section C.5 introduces a change of variables for the integral equation. The operation 
provides a form for the integrals allowing them to be made independent of the com- 
pressibility factor, p. 

Section C.6 derives the transformation from the compressibility coordinate system to the 
local coordinate system for each panel. Applying this transformation to the integral 
equations provides the panel integrals which are evaluated in app. D. 

C.l PRANDTL-GLAUERT TRANSFORMATION 

The Prandtl-Glauert transformation is used in the following to relate solutions to Helm- 
holtz’s equation with solutions to the unsteady flow equation, eq. (30). Using 
indicial notation and summation convention (ref. 5, sec. 7), the transformation is 
expressed as follows: 

x* = Ejixj (C.l) 


t This distinction in measure and its consequence, viz., the introduction of two different 
unit vectors which are both normal to the same surface, tends to cause confusion in 
following the development of the panel method; it should be carefully noted. 
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where the compressibility coordinates are denoted as 


1 2 
X =x,-x^ = y, 


x^ = z 


and 


® 1 ~ ^ ’ ^2^ ~ ^3^ ~ ~ 0 otherwise; 

and the Prandtl-GIauert transformation is introduced choosing the x * system as the 
reference frame (ref. 5, pg 8). The metric tensor (ref. 5, sec. 29) of the compressibility 
system, therefore, is given by 

9x^ 9x^ 

8ij= — : r 

^ 9x* 9xJ 

where, by eq. (C,3), 

^1 ] ” ^ ’ ®22~ S33 “ gjj = 0 otherwise; 


while, its contravariant form (ref. 5, sec. 30), viz.. 



g 


where G^J are the cofactors of gy and g is the determinant of gy, 



has the components 


gl 1 _ 1 ^ g22 _ j/^2 g33 _ 1/^2^ gij _ Q otherwise. 


Helmholtz’s equation (ref. 3, pg. 239), viz., 

V 2 </,* = 0 , 

is expressed in terms of the contravariant metric tensor, eq. (C.6), by introducing 
the following identity (ref. 5, eq (90.7)): 

V2 0*=g'j0*,y 


(C.2) 


(C.3) 


(C.4) 


(C.5) 


(C.6) 


(C.7) 


(C.8) 


(C.9) 


(C.IO) 
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where (ref. 5, see, 33) the comma denotes covariant differentiation. For the case of 
linear transformations, such as the Prandtl-Glauert transformation, covariant differenti- 
ation reduces to partial differentiation; hence, . 



9^0* 
9x^ 9xj 


(C.ll) 


Substituting eq, (C.l 1) into eq. (C.IO) and substituting the result into eq. (C.9), Helmholtz’s 
equation becomes 


9^ 0=* 


9x^ 9xJ 


- + J2"^0* = 0; 


(C.l 2) 


and, introducing eq. (C.2) and eq. (C.8), this form of Helmholtz’s equation is expressed as 


Helmholtz’s equation, as given by eq. (C. 13), when multiplied by is identical with the 
unsteady flow equation, eq. (30). As a consequence of this identity, solutions to Helmholtz’s 
equation expressed in the compressibility coordinate system (regarding the scaled coordi- 
nates, x\ as the reference frame) are solutions to the unsteady flow equation. 


In the scaled coordinates, x\ Helmholtz’s equation appears as follows: 


9^ 0* 9 — 

+ 0* = 0; 

9x*^ 9x^ 


(C.14) 


or, letting 


Helmholtz’s equation becomes 


— - 1 - -2 - -3 

x = x ,y = x"^j z = X , 


0 --* + 0 --* + 0 --* + n^(j)* = 0 


(C.15) 
(C.l 6) 


If 


0*(x,y,z) = 0*(x,|3y,i3z) = 0*(x,y,z) 


(C.l 7) 
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where 0 is a solution to eq. (C. 1 4), then, on substituting eq. (C. 1 7) into eq. (C. 14), the chain 
rule for differentiation leads to 


3^0* 9x* 
3x^ 3xj 9x^ 


9xJ 9 
• — - + 0* ^ 0 . 


(C.18) 


Introducing the inverse Prandtl-Glauert transformation, viz., 

x^ = x-j 


(C.19) 


where 

S j 1 , ~ ~ 0 otherwise, (C.20) 

eq. (C.18) becomes eq. (C. 13); hence, if 0* is a solution to Helmholtz’s equation expressed 
in the scaled coordinates regarded as the reference frame, the change of variables shown 
by eq. (C. 17) provides a solution to the unsteady flow equation. 


C.2 HELMHOLTZ’S THEOREMS 

Solutions to Helmholtz’s equation are provided by the theorems of Helmholtz which 
appear in the following. Since the preceding has shown that solutions to Helmholtz’s 
equation can be related to solutions of the unsteady flow equation by the Prandtl-Glauert 
transformation, the following theorems provide a means for constructing solutions to 
the unsteady flow problem. 

Helmholtz’s theorems (ref. 3, pg. 240) are based on a division of all space into two volumes 
V and V' by the closed surface 2 enclosing V; thus, P, an arbitrary point in space, is 
either interior to 2 (i.e., PeV) or exterior to 2 (i.e., PeV') or on 2 (i.e., Pe2). Letting 


47tR 


(C.21) 


where R is the distance from P to a point Q on 2 and noting that eq. (C.21) satisfies 
eq. (C.12), then, if 0* is a second solution, the first theorem of Helmholtz states that 


'//'(IS (‘ 5 ’') 

2 

0*(P) = OifPeV' 

0*(P) = oo if Pe2 


(C.22) 


125 



where ^ (Q’*') denotes the limiting value of 0* as the point of evaluation approaches Q 
from V, n' is the normal coordinate of S and is positive in V, and the prime denotes 
operations involving the coordinates of Q. The second theorem states that 

0* (P) = -JJ (Q") 0* - 0* (Q') ] ds if PeV' 

S 

= 0 if PeV (C.23) 

= oo if PeS 

* - . * 

where 0 (Q~) denotes the limiting value of 0 as the point of evaluation approaches 
Q fromV'. 


Combining eqs. (C.22) and (C.23) leads to a third theorem: 



0 *- 110 *] 


30*\ , . 

i?r ‘ 


if PeV' UV 


where, for example, 


(C.24) 


110*1 = 0*(q+)- 0* (q-) (C.25) 

because 0 and 90 / dn' are not finite at the surface. The elementary source, 

0*, leads to a discontinous normal derivative for the potential; while the elementary doublet, 
90*/9n', leads to a discontinuous potential. Noting that 0* and 90*/9n' represent, respectively, 
the complex amplitudes of the potentials induced by unit sources and unit doublets distributed 
on 2, one may define 


and 


a* (0) = 



(C.26) 


M* (Q) = tt0*l 


(C.27) 


so that 


0* (P) = 0g* (P) + 0j)* (P) 


(C.28) 


126 



where 


<j>* (P)^ JJa* xIj* ds' 


and 


0D* = " 


-If 


Further, the complex amplitude of induced velocity can be expressed as 


V* (P) = V0* = T * (P) + Vtj* (P) 


where 


and 


because 


Vg* (P) = -JJ J* 
S 


v' \p* ds' 




ds' 


v\p* = - v' ||/* 


where the prime denotes differentiation with respect to the coordinates of the point Q, 
while the unprimed gradient operator denotes differentiation with respect to the 
coordinates of the point P. 

Eqs. (C.28) through (C.33) describe the complex amplitudes of the potential and 
velocity induced at the field point P by harmonically fluxuating sources and doublets 
distributed on the closed surface 2 . These induced quantities satisfy the form of 
Helmholtz’s equation shown by eq (C. 1 2); and, eqs. (C.28) through (C.33), like 
eq. (C.l 2), are valid in terms of any coordinate system which is the result of an 
admissible transformation (ref. 5, sec. 19 and 20) from the reference frame. The 
reference frame has been chosen, in sec. C.l, to be the scaled coordinates, x*. 


(C.29) 


(C.30) 


(C.31) 


(C.32) 


(C.33) 


(C.34) 
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C.3 CONORMAL VECTOR 

As noted in the preceding section, the integral equation, i.e., eq. (C.24), which provides 
a means for evaluating the unsteady flow problem, is expressed in terms of geometric 
quantities whose measure is that of the scaled coordinates, This use of measure is 
illustrated by writing 

90 * /s 

— ds = n • v0* ds 

9n 

where the magnitude of n is unity in the measure of the scaled coordinates and ds 
is the differential surface element in the scaled coordinates. The objective of this 
section is to derive the following identity : 

n ds = lij, • v0* d S 

where rTg is the conormal vector with components 
where the components of N satisfy 

^(N,)^ + (N2f + (N3)^=I 

in the measure of the compressibility coordinates, x^, and dS is the differential surface 
element in the compressibility coordinates. The vectors, n and N, are both normal 
to the same surface ; hence, when they are evaluated at the same surface point, 

n = N/N 

where N is the magnitude of N in the measure of the scaled coordinates, 

Further, the two differential surface elements are shown to be related as 

dS = N ds . 


As a preliminary to the derivation, consider the reciprocal base systems (ref. 5, 
sec. 45). In the reference frame, again, taken to be the scaled coordinates, x\ 
position is described as follows: 

r = x^ bj 


(C.35) 


(C.36) 


(C.37) 


(C.38) 


(C.39) 


(C.40) 


128 


where (ref. 5, eq. (45.2)) the vectors, bj, are the unit base vectors of the scaled coordinates 
such that 

(C.41) 

The covariant base vectors of the compressibility coordinates (ref. 5, eq. (45.6)) are 
defined as 


9r 

SL: “ ' I 1 

1 3x* 


(C.42) 


and, since 


aj-aj-gij , (C.43) 

these are not unit vectors in the measure of the scaled coordinates. 

Contravariant base vectors are introduced by eq. (45.8), ref. 5, and these can be written 
in the follovwng convenient form: 


aWyk=ljXaj, (C.44) 

where (ref. 5, pg. 138) 

^ijk=V?eijk (C. 45 ) 

and ejjk is the permutation symbol defined in sec. 40 of ref. 5. Similarly, from 
eq. (4^9), ref. 8, 


where (ref. 5, pg. 1 38) 




e ijk = _L gijk 


(C.46) 


(C.47) 


Returning to the derivation, as noted by eq. (90.5) of ref. 5, 


d(p* 

— ds = n 
3n 


i 


30* 

— r ds 
3x* 


(C.48) 
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where denotes the contravariant components of the unit normal vector, n. An expression 
for evaluating n^ is obtained from the description of the surface, 2, in Gaussian form 
(ref. 5, sec. 52), viz.. 


x^ = x^(u^,u^) 

where u®' (for a =1,2) are parameters which, when varied independently, describe 

curves on the surface which are surface coordinates. The base vectors of the surface 
coordinates (ref. 5, eq. (54.6)) are given by 

a"? 




(C.49) 


(C.50) 


where r" denotes position on the surface, i.e., 

"?= X^ (u^,u^)bj = ^u^,u^)"aj 

where the following has been used : 


(C.51) 


so that, by eq. (C.40) and (C.20) , 


and 


^ dr 9xJ 


bi = 


9xJ 9x* 




fij X* (ul,u^) =x-i (u^u^) . 


Since a’l anda ‘2 are vectors tangent to the surface along the surface coordinate curves, 
the unit normal vector is given by (ref. 5, eq. (65.2)) 


A 

n = 


^1 X ^2 
^1 ^2 


(C.52) 


(C.53) 


(C.54) 
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The components of the unit normal vector are derived using eq. C.50 and C.5 1 to write 
the surface coordinate base vectors as 



9x^ 

9u“ 



Their cross product, using eq. (C.44), is then expressed as 


3] Xa2 


dxj 9x^ ^ X 9xj 


3u^ 3u^ 


hence, defining 


/a= ajx'a2 , 


the covariant components of the unit normal are found to be 


= 9x' 3x^^ 


while the contravariant components are 


= gy nj 


(ref, 5, sec. 45). 


As shown by sec. 54, ref. 5, the differential surface element is given by 

ds=|'?jX‘a2 du* du^ ; 

and, introducing eq. (C.5 7), 


ds = du^ du^. 


Combining eq. (C.58) and (C.61), the covariant components, viz.. 


3x-i 3x^ 1 

' -1 Hn ^ ( 


"i^s=^iik“^ du\du“ , 
au’au2 


(C.55) 


(C.5 6) 


(C.57) 


(C.58) 


(C.59) 


(C.60) 


(C.61) 


(C.62) 
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are seen to be the components of the vector formed by the cross product of the following 
two tangent vectors: 

dr| = aj du^ = dxj * aj 


and 


dT2 = dx2^Hi 


Consider the following vector product: 

N dS = dr j X dr2 


in the measure of the compressibility coordinates; thus, 

N dS = ejji^ dxjj dx2'' . 

Introducing eq. (C.62) and recalling eq. (C.45), eq. (C.65) leads to 

N dS = Vg nj dS a^ . 

Taking the norm of this result in the measure of the scaled coordinates, wherein 

^nj'ai^ • (nj'ai^ = nj g^ nj = 1 

and 

(Njti) • (Njti) = Nj g« Nj = (N)2, 

the two differential surface elements and the covariant components of the two normal 
vectors are found to be related as follows: 


and 


ds = Vg N dS 


n-Nj/N, 


Eqs. (C.69) and (C.70) verify the expressions shown as eqs. (C.38) and (C.39). 


(C.63) 

(C.64) 

(C.65) 

(C.66) 

(C.67) 

(C.68) 

(C.69) 

(C.70) 
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Recalling eq. (C.59) and introducing eq. (C. 70), 

n^ = .g^ Nj/N ; 


thus, since (ref, 5, eq. (92.5)) 


d(p* : d<b* 

— ds = n^ ds , 

3x' 


eqs. (C.69) and (C.7 1 ) lead to 


which can be written as 


where 


90* 

9n 


ds = VI 


g‘JNj 


90'‘ 

9x* 


dS 


90* 

9n 


1 i 

ds = n^,i— dS = 
9x* 


30* 

9nc 


dS 


n^'=VI 

and, on introducing eq. (C.8), (C.5), and (C.7), 

= ^j^^2 = N2,n/ = N3. 

From eq. (C.76) the vector n^. is seen to be the vector usually called the conormal to 
the surface. 


The integral equation shown as eq. (38) follows directly from eq. (C.25) by making 
the substitution shown as eq. (C.74); thus. 


where 


0* (P) = 


//(• 


0* - m'* 





(C.71) 


(C.72) 


(C.73) 


(C.74) 


(C.75) 


(C.76) 


(C.77) 


(C.78) 
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C. 4 ELIMINATION OF PANEL EDGE LINE VORTICES 

The integrals shown by eqs. (C.29) through (C,34) are expressed by the sum of the 
panel integrals shown by eq. (37) through (60). In the case of eq. (60), viz,, 

Vd* (P) = - ^7^, (v' <//*) dS', 
a decomposition is introduced such that 


(C.79) 


where the term represents the velocity induced by a line vortex lying on the perimeter 
of the panel, viz., 9Sg. When panels are joined together to form a continuous surface of 
integration, the line vortices at the edges of panels ^abutting one another are equal in 
strength and opposite in sign; hence, the velocity v^ induced by one panel edge is canceled 
by the velocity induced by the other, abutting, panel edge. As a result, the line vortex 
induced velocities need not be computed and the doublet panel induced velocities (i.e., 
eq. (60)) are computed from the regular part of eq. (C.79), viz. vj . 


The decomposition of eq. (60), which leads to eq. (C.79), is formed from a corollary to 
Stake’s theorem. Stake’s theorem is expressed by eq. (92.15) of ref. 5 which, on inter- 
change of dummy indices j and k, is written as follows: 

Fkj "i ^ '''= 

2 92 

where 92 denotes the perimeter of the surface, 2 ; 

x^ = x^ (c) 

represents a parametric description of 92; and Fj^, represents the covariant com- 
ponents of an arbitrary vector. The corollary is obtained by letting 


whereby, eq. (C.56) becomes 


dx‘< 


=*km s'”" («♦ *r„) . "i <1» = ff =ekm ^ 


(C.81) 


92 
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where use has been made of Ricci’s theorems (ref. 5 , sec. 35), viz., 

gj.,k=0 and g'i,^ = 0, 

whence 

(^8km),j ®Ckm ~ 2Vg ®£km ~ ^ • 

Expanding (C.81) leads to 

Ata K ^ ^’"""5 < "i 

as 2 

s 

The vector form of the integrand of the second integral is deduced using the reciprocal 
base system shown by eqs. (C.44) and (C.46). The result is a triple vector product which 
is then revised using the following vector product identity: 

("a X 1b) X "c = b (^ • "c) -"a (b • ^). 

These operations appear in the following sequence : 

egkm g™" Kj eUk (tk X A,) nj 

= - p* (“a* xtO X “aj^^ g*^" A Uj 
= - ju* (n xv) X V i//* 

= - M* (n ’ v)vi//* + ;u*(v • V 0*) 

= - M* g‘J nj i/',gj*'a^ + M* ng gh (i/,jj*'a® . 

Also, from eq. (C.12), 

g’-’ = -n~ \p*; 


(C.82) 


(G.83) 


(C.84) 
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and, because the first term is identical with integrand of eq. (60), viz., 

ju* n- = M* — (vi//*) , 


eq. (C.79) becomes 


-1.* ->.* 
VD =Vd +v^ 


where 




= J'J^ M*,j ^*,^1 nj ds'a*^ + ju* ds'^a® 


= J^y*(n X vju*) X v' il/* ds' + ^2^ n ix* \p* ds' 


% =- / 


p „mn . * 

^«km ^ 


dx r 

— dc a^ = / /x*drxv 

dc J 


m 


9S, 


92. 


where 


and 


>i ~ ’i' > 'Z'* . 


dx ^ 

dr = - — dc au . 
dc ^ 

C.5 CHANGE OF VARIABLE 

In section C.l it was noted, eq. (C.18), that solutions to the flow equation could be 
related to solutions to Helmholtz’s equation through a change of variables: thus, if 
0* (x,y,z) is a solution to the flow equation, then 

0* (x,y,z) = 0* (x, y/|3, z/|3) 

is a solution to Helmholtz’s equation expressed in the scaled coordinates, viz., 

^x*x ^y*y ^zz 0* = 0 • 


(C.85) 


(C.86) 


(C.87) 


(C.l 7) 
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The advantage accrued through this change of variables, when the solution to the flow 
equation is the integral equation shown by eq. (C.77), is that the integrals do not contain 
the compressibility factor, j8, except that the compressibility factor appears in the 
Helmholtz’s coefficient, viz., 

n = <3 Mq ^ (3Q) 

* 

which appears in the exponential function contained in \jj' . 


Recalling eq. (C.77). the panel induced complex potential amplitudes at the field point, 
P, are given by 


and 



(C.88) 


(C.89) 


Under the change of variables eqs. (C.74), (C.72) and (C.59) yield the expression 


9i//* 

9"c 


dS = njgy — i dS 

^ OX 


which becomes 


^c' 


9i//* 

dS = n:— dS 
>9x ‘ 


(C.90) 


where 


nj nj 


9i//* j 


and 


hence, the doublet panel integral becomes 


(C.91) 


137 



where 


when 



L =- L 

9n ^ 9x^ 


(C.92) 


Recalling eq. (C.69), the source panel integral becomes 


0* (P) = J (P) 

where 

1 __ 

VgN 


(C.93) 


and 


0g* (P) ~ • 

Having the potentials expressed in the scaled coordinates, a convenient means for 
computing the velocity components in the compressibility coordinates is obtained 
by the following operations: 


* _ ^ £ j 

9x‘ 9x’ 9x* * 9x-> 


(C.94) 


(C.95) 


Noting, again, that 


the source panel induced velocity components are given by 


Vgi (P) = JEJvsj (P) 


where 


'si • 

•1*' 


The doublet panel induced velocity components are given by eq. (C.85) whereby 


(P) = 


n-ds'.«2 


9xJ 9x 


JJ'^9. iu* i//* 


i//* ds' 


so that 


VDi(P) = Ei«voj(P). 

The influence of a panel edge line vortex, eq. (C.86), becomes 

_» r dx^ 

VvB (P) = 7 M* — dc 

2, 


9^* 
9x'« ■ 


C.6 PANEL, LOCAL COORDINATE SYSTEMS 

Eqs. (C.91) through (C.lOO) of section C.5 contain the integrals which must be evaluated 
to determine the panel influence at field points. This section introduces local coordinate 
systems in which these integrals assume forms which are particularly convenient for that 
evaluation. The coordinate lines of a local system are denoted as rj, f and the system is 
defined such that the t? coordinate plane is parallel with the panel plane and such that 
the f coordinate is along (and in the direction oO the panel normal, n. 

Letting the coordinate lines be denoted as 

tbe local system is expressed by a transformation of the compressibility coordinates 
as follows: 

^i = Ajixj 


(C.96) 


(C.97) 


(C.98) 


(C.99) 


(C.lOO) 


(C.lOl) 


(C.102) 
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with the inverse; 


This is a linear transformation since it is a combination of the Prandt-Glauert transforma- 
tion, viz., 

jfk _ ^ 

and an orthogonal transformation, 
such that 

Ad = tl ^ E-^ 

Aj tg iij 

where, because of the orthogonal nature of the transformation, 
and 


The orthogonal transformation is derived by noting that, if oij are unit base vectors 
for the local system, then they are related to the base vectors of the scaled system by 
the transformation rule for covariant tensors (ref. 5, sec. 24), viz.. 


A 



3^^ 


— ~ =TJb= 

3 x-> 


The unit base vectors of the local system are chosen as follows: 


«2 



A A A 

a I = 0!2 X n 



(C.103) 


(C.l) 


(C.104) 


(C.l 05) 


(C.l 06) 


(C.l 07) 


(C.l 08) 


(C.109) 
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where n is the unit vector normal to the panel plane defined by eq. (C.54) and c is a unit 
vector in the direction of the freestream flow, viz., in the direction of x and x. Because of 
the form of the Prandtl-Glauert transformation, the vector c is a unit vector in the measure 
of both the scaled system and the compressibility system, and its covariant and contra- 
variant components are identical to one another. Using these characteristics of the vector c 
along with eq. (C.38), the unit base vectors of the local system are seen to be given by 


and 


«2 


A 



Nxc 


X cl 


In 


1 A ^ 

— X N 
N 2 


«3 = 



(C.llO) 


where it should be noted that these vectors can be formed directly from vectors expressed 
in the measure of the compressibility system. 


The expansion required by eq. (C.108) is obtained by writing the cross products using the 
reciprocal base systems (viz. eqs. (C.44) and (C.46)) ; thus. 


N X c = e^-i^ Nj Cj^'a^ . 

Introducing 

v^ = ey^NjCj, 

and 


(G.lll) 


(C.112) 


a = N X c 




Rii v-* 


(C.113) 


because the local system has the measure of the scaled coordinates, wherein 



the unit base vectors of the local system are found to be as follows: 
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(C.114) 


a: 

^ a . 


and 


A 1 \ ^ 

a.- — U^ a: 
1 aN ^ 


=rr N,' a-^ =77 N; g-^^’a: 


N J N J 


where 

ui = eij*^ Vj Nj^; 

because, in the measure of the compressibility system, 

vi = Vj. 

Since, again, by the law of transformation for covariant tensors. 


37 37 3x-i i A 

a =—=_ — = eJ b:, 

3x* 3x‘3x‘ ^ 


the local system base vectors become as follows: 




-V* EJ b: 
a ^ J 


a, = — uiEjJ bj 
‘ aN ‘ J 


a- 


— Ni E-J b- 


The orthogonal transformation shown by eq. (C.106) is deduced from eqs. (C.108) and 
(C. 1 1 8) with the result that 


1 


xJ - U> EjJ ^ V- vi E J ^2 + i N. gki E j ^3 
aN a ‘ N • ‘ 


(C.115) 


(C.116) 


(C.117) 


(C.118) 


(C. 1 1 9) 
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Similarly, the transformation, eq. (C.103), which is the inverse of eq. (C.102), is found 
from eq. (C.l 14) as 

xi = ~ +TrNigj^^^. 

ckN ^ OL N J 

Using the characteristics of the orthogonal transformation shown by (C. 1 07), the trans- 
formation shown by eq. (C. 104) is found as 


and 


= ^uiEJ’xj, 

aN ^ 

^2 ^ IvifiJ xJ , 

O! V 

j3 


Finally, introducing eq. (C.l) into eq. (C.l 21) and noting eq. (C.4), viz.. 


and noting from eq. (C.6) that 


<j.. = E.^ E-^ 

1 ’ 


g« gi“ = Si*^ 


the transformation shown by eq. (C.102) is found to be 


gjj xJ , 
OL 


and 


o 1 : 

— N; XJ . 
N J 


♦ 




(C.l 20) 


(C.121) 


(C.l 22) 
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Matrix notation provides a useful form for. presen ting the operations involved in computing 
the components of Aj^ and its inverse, viz., 




(C.123) 


and 


where 





(C.124) 


{uH = components of (N X c) X N , 

{v*} = components of (N X c) and (C.125) 


Nj| = components of N. 


Expanding eq. (C.95) and noting eq. (C.105), the velocity components are found to 
transform as follows: 


V;* = EJ^ 

ax^ 


■ 3xi ‘ ' 


k^* 


A,k 


where 0 denotes the complex amplitude of the modified potential expressed in the 
local system and 



(C.126) 


(C.127) 


denotes the components of the velocity in the local system. Using eq. (C. 1 26), eqs. (C.97) through 
(C. 100), describing the panel influences, become as follows: 

(1) Source panel induced potential: 

(P) = J0 g* (P) (C.128) 

where 


(P) 


=-//. 


* 0* d^^di?^. 


(C.129) 
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(C.130) 


(2) Doublet panel induced potential: 


</>d(P) = ?d(P) 


where 


dr dr? 


(3) Source panel induced velocity components 


^ . j 

^Si = V^Sj 


where 


V 

Sj 


// 


- , , 

d* — dr dr?' . 

ar' 


(4) Doublet panel induced velocity components 


VDi = 


where (for Greek indices ranging over 1 ,2) 


because the components of the normal are given by 

nj = 0, n2 = 0, 113 = 1 


and 


while 


^8km = - ®S/3m = ' ^5/3 


7d3 (F) = , „2 J. 


ar ar 


(C.131) 


(C.132) 


(C.133) 


(C.134) 


(C.135) 


(C.136) 
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because 


„3jk g = gSaiS = paj3 „ „ 

® ®3km ® ®3j38 ® ®|36 


(5) Doublet panel edge induced velocity components: 


V*.= A ^ v*o 


(C.137) 


where (for Greek indices ranging over 1 ,2) 


\5 (P) 




■'5a' 


d^“ 9i//* 

dc —rr 

dc ar 


as. 


(C.138) 


while 


d^“ . a^* 


vv*3m = - 

as„ V 


(C.139) 
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APPENDIX D 

EVALUATION OF AERODYNAMIC INFLUENCE COEFFICIENTS 


This appendix describes the theory underlying the procedure used to evaluate the panel 
influences shown by eqs. (C.l 28) through (C. 136) of app. C. 

Panel singularity distribution functions are introduced in sec. D. 1 , These functions reduce 
the panel influences of app. C to influence coefficients relating the induced flow potential 
and velocity at a field point, P, to the coefficients of a T ay lor series expansion of the 
distribution function, cf. eq. (49) and (52). The end result is a set of integral equations for 
the panel influences which are in terms of panel integrals having integrands of nine basic 
types. The remainder of the appendix is directed toward deriving evaluation procedures 
for those nine integrals. 

Section D.2 introduces two coordinate systems: (1) a panel edge coordinate system describing 
position in the panel plane in terms of distance along and normal to a panel edge and (2) a polar 
coordinate system. In the following sections of the appendix the panel integrals are expressed 
in polar coordinates. The panel edge coordinates are then introduced as a change of variables. 
The result is a convenient expression for the integral in terms of the distance along the panel 
edge. 

Section D.3 contains a reduction of the integrals involved in the panel influence coeffi- 
cients. The reduction is based, in part, on a separation of the unsteady kernel function into 
a steady component and an unsteady component. The unsteady kernel function, viz., 



47tR 

is expressed as 

\jj* = + *\jj 

where the steady component is given by 

1 


and the unsteady component is given by 

47t ' R ■ 

The steady component is evaluated in closed form and the unsteady component is further 
reduced to a dependence on two fundamental integrals which cannot be evaluated 
directly in closed form. 


The integrands of the two fundamental integrals are approximated in sec. D.4 by polynomials 
containing terms which are closed form integrable. Sec. D.5 describes the polynomial 
estimation procedure for the approximation introduced in sec. D.4 and sec. D.6 contains an 
error analysis to the approximation. 

Finally, sec. D.7 contains a derivation of the influence coefficient evaluation procedure used 
in the far field. The fundamental idea underlying the far field evaluation is that of approxi- 
mating the kernel function (or its gradient) with a truncated Taylor series expansion which 
is a valid approximation when the field point is sufficiently far from the influencing panel. 

The section presents the criteria used to establish validity. 

D.l PANEL SINGULARITY STRENGTH DISTRIBUTION FUNCTIONS 

For the purpose of introducing the singularity distributions, o*(Q) and m*(Q)> the panel 
influence integrals are expressed in terms of the following representation of the field 
point, P, and the surface point, Q: P' is the point representing the projection of the 
field point onto the panel plane and h is a vector along the panel normal from P to P', 
figure 50. A local ori^n , Q^, with coordinates rjg, fQ_is introduced and P' and 
Q have the positions 8 and "r relative to Qq. The vector R is then given by 

R = p + h. (D.l) 


P 



Figure 50. — Representation of Field and Surface Points 
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Indical notation and summation convention (ref. 5, sec. 7) are used in the following; 
hence, the following notation is introduced for the coordinates of points: 


.1 


= v' : coordinates of Q 
^o “ ^o > ^o “ ’^o • coordinates of Qo 

1 9 

z = ^ , = T? ; coordinates of P' 


The components of the vectors in the panel plane (viz.,"?, and S) are given by 

where, and in the following, Greek indices have the range a = 1,2. 

The panel singularity strength distribution functions, c.f. eqs. (49) through (52), are 
expressed in terms of the above notation as follows* : 


a* (Q) = a* + a* ^ r« 


and 


“• «) = f‘i + , op r“ / 4 / r\ 


For the purpose of evaluating the^integrals of eq. (C. 1 28) through (C. 1 36), however, it is 
convenient to expand a and fx in a Taylor series about the point P' (viz., the 
projection of the field point onto the panel plane); thus, 

a*(Q) = a^ + a^p« 

and 

«• (0) =■ i + i / / 


* Terms which are of third degree in the spatial coordinates are included in the 
doublet distribution function for evaluating the quasi-near-field approximation 
described in sec. D. 7. 


(D.2) 


(D.3) 


(D.4) 


(D.5) 
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where 


III 

o 


III 

O 

^o,a 

III 

o 

=l 

M$+PS,a8“ + M*,a^ 
Po,a + <a^S^ + ^^^ 


^O!0 “ 

Po,o:i3 

3 

III 

^o,<x^y ■ 


Introducing eq. (D.5), eqs. (C.128) through (C.136) are written as follows: 

0^ + ix^oj (p“) + ^ CO (p“ p^) + (p“ P^)j 

hw^O) + a^W^(p«)] 

Vg*3=[^aoco(l) + ffaw(p®') 

^da " (J-a(^0) + Mq,|3 co (p^) + co (p^ p'^) 

L ^ 


where 


4 = - [fc, W,, (1) + M„ J (/) W„ (pO ^'1')] 

+ n2 Lj, * (l) + n^ J, (<’“)+|^“a(3 * 


0 (f)=//0* f d^' d??' 
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(D.14) 


w„(f) s/Jg-ll fdl'd.,' 


The problem of evaluating the panel influences reduces to the evaluation of the 
following integrals: 

o> (1 ) , CO (p“), 00 (p“ p^) , w (p“ p^ pT) 

Wq, (1) , Wq, (p^), Wq, (p^ p^) . 

D.2 PANEL EDGE COORDINATE SYSTEMS 

The corners of a panel are enumerated by the subscript e in a counter clockwise 
direction. Their positions are denoted The eth edge of a panel is described by 
the vector 

"^e ~ ^e+.l ” ’ 


and 


^e = fe/|f,| 


is a unit vector tangent to the eth edge. In indicial form 


where denotes the unit base vectors of the local system, cf. eq. (C.109). A normalized 
edge normal vector is then introduced as 


A 



®ccj3 ^e 


a 




such that 
and 




X L = 


a 


' 0!|3 




- n, 


2, 1 


t.* = 1 


(D.15) 


(D.16) 

(D.17) 

(D.18) 

(D.19) 


151 


The orthogonal unit vectors n^ and U are used to write 


whence. 


P^agng + Vte 


^ A 

dp = tp dc 


(D.20) 


(D.21) 


where dc is the differential arc length along an edge; also, 

dc = dV. 


(D.22) 


Recalling eq. (D.l) and introducing eq. (D.20), 


R = ag ng + vtg + h ; 


(D.23) 


whence, 


^ 2 = 2 ^^2 


Introducing polar coordinates, viz.. 


1 9 

p = p cos (/» and p^ = p sin 


(D.24) 


where 


I PI- 


then 


p X dp = p X 


P / ^ 

— dp + p d 0 ( - sin (i) cei+ cos 

P ' 


<A«2) 


Letting u = p d<j> sin cx I + cos 0 a- 2 ^ 


p X eQ,g = p^f^ - p2 i; ^ = p^ d 0 ^cos^ 0 + sin^ p2 d 0 ; (D.25) 
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thus, 


-A. 

p X dp = p- d(j) 


Alternately, 


P X dp = ("ag Hg + vtg) X tg dv = ag dv ; (D.26) 

hence, combining eqs. (D. 25) and (D. 26), 

ag dv 

d0 = - (D.27) 

3gZ ^ yZ 

provides a relationship between the differential of the phase angle, 0, and the edge 
coordinates. 

D.3 REDUCTION OF THE PANEL INTEGRALS TO FUNDAMENTAL n^TEGRALS 

Proceeding to a consideration of the integrals to be evaluated, consider fii^t the relationship 
which can be established between co(p“0 and Wgj,(f). 

Since 


d\jj* _ ^ 

9R 


(L 

\R 9R / 



) 


fd^' dr?' 


while, on the other hand, since 


it follows that 



Consequently, the two integrals are related as 

, w(p«f) = -hW^(f); 
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and, by induction, 


co(p«) = -hW„(l) 
w(p“p^) =-hW^(p/^) 
co(p“p^ pT") = -hWo, (p^ p'y) 


The remaining integrals (viz., \p, co, W^,) are evaluated in terms of the cylindrical coordi- 
nates; whence. 


d^ dp = 


9 , P) 

d (p ,(p) 


dp d0 = p dp d0 . 


In general, these integrals are not convergent; rather, their correct evaluation must be 
performed using the theory of finite parts. 


Consider the integral performed for each edge, viz.. 



Pe (0) 


I 


gp dp 


where (figure 51) and denote the coordinates of the end points of the edge, 
while pg(0) is simply the value of p expressed as a function of 0 along the eth edge. 



Figure 51. — Region of Integration Corresponding to a Pane! Edge 


(D.28) 


(D.29) 
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Using these concepts, 


Defining 


and using the notation 


it follows that 


<Pe PqW 

i//(l) = S J" d(j) j i//*/odp 
<Ae" ° 


R 


X(R)= j i// R' dR' 
o 


Re (</>) =p (0) , 


0p 




^ ^ ^ e / d0 [X(R)-X(,|h|)] =2 J d0 X (R) - 27 t Cq X 


(|h|) 


<t>e 


(pa 


where denotes the winding number, if the origin of the p^ coordinates is interior to the 
panel, 2g. The winding number is defined by the following relation; 


f d0 = 27 t Co = 

* 


o (the boundary of 2g does not 
encircle the origin) 

2ir (the boundary of does 
encircle the origin 


Consider, now, the integral 


<1' (p“ f) = 0* f d^' dp' (r rp*) (0 dt' dp' 


(D.30) 


(D.31) 


(D.32) 


(D.33) 


(D.34) 
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Integrating by parts using Gauss’ theorem, 


^P(p^f)= fn„Xfdc'-/y*X— d^dTj'; 

J JJ 


thus, by induction, 


li/ (p“) =J X dc' 

a 

= Jpl^n^Xdd -JJxd^' dr}' 
a 

^ (p“ P^ P'^) = fp^ n^ X dc' - X d|' dr?' - ^ dr?' 


where 


= I if a = ^ 
= 0 if . 


The integrals 


JJxd^'dr}' and //^ “ X dt dr? 


(X Y At' A^' 


may be treated in a fashion identical to that given to i//(l) and \p(p^). 
Defining 

R 


9 = J X R' dR' 
o 

<t>e^ 


and 


X(l)=J^XdrdT?' = 2 J 0d0-27rCo0 
X (p«) ^jjp X dr dr?' - y* 0 dc'. 


(Ihl) 


(D.35) 

(D.36) 

(D.37) 


(D.38) 

(D.39) 

(D.40) 
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Turning now to the treatment of co(l) and Wq,( 0 and recalling that 


9R h Cfb\b* 

5p=--. o.(i)=j|/^ df'd,' 




_ _ 1 3;//* 

=-h/i-- drd.'=-h 




Pq (<l>) 


?/ i 




P dp 


But, 


hence. 


p dp = R dR ; 




C0(l) = -h| ^ y \i/*d0-27rCoi//*(lhl) 


Next, 




0* f ds' - 



dr dT?' ; 


consequently, by induction. 



where the integrals on the right are given by eq. (D.28) 


(D.41) 

(D.42) 

(D.43) 

(D.44) 
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All of the required integral evaluations are now seen to be in terms of the following 
two operators; 

the panel operator: , 


J(f)=E f fd(/.-27rCof(|h|) 


(D.45) 

6 

^e' 

the edge operators: 



le (g) = y* g ds' 
E, 


(D.46) 

\H1) = J(X) 


(D.47) 

^ (P“) =Ele(X)(ne)„ 
e 


(D.48) 

</' (p" P^) =E le (P^ X) (Hg) ^ - J ( 0 ) 

e 


(D.49) 

*(AV) = £ le(A'‘) (n,)„ - 6p„ I, (9) (0) ,n^)^ 

e e e 

(D.50) 

w(l) = -hJ(if*) 


(D.51) 

W<,0)= Y, ‘eWtCe);, 

e 


(D.52) 

W>(i)= £ I^(A*)(n,)„ -*(J) 

e 


(D.53) 

W« (p V)= x: ^e (p *) (ne)a ~ (P^) ' 

e 

(p^) 

(D.54) 

A careful examination of the integrals appearing in eqs. (D.47) through (D.54) reveals that the 
only integrals which need be evaluated are the following: 


Ie(p"p^i//*) 


J(X) Ig(X) Ie(-°"X) 

Ie(pVx) 

(D.55) 


Before attacking this list of integrals, consider explicit expressions for the functions X 
and d, viz., 


IX IX 

X = -47r J xp*R'dR' = J e'‘^*^dR' = - 


0= / XR'dR' = 


J R' -l)dR' 


= -i- {(-iJ2R) [e-iJ^R- 1 -(-iflR)] 

(iJ2)2 

- - 1 - (-i«R) - (-iJ2R)2 / 2 ]} . 


Introducing the abbreviated notation 


X ■ -iJ2R , 


-4ir\p* = / R ; 


X = R 


'e^-i 


0 = R^ {(e^ - 1 - x) / x^ - (e^ - 1 - x - x^ / 2) / x^| 


From these expressions 0 can be seen to be related to \p* and X as 

n-d =4nR^\P* + Rx / 2 + X 

= 47rR2)//*-iJ2R2/2 + X. 
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Returning to the problem of evaluating eq, (D.55), recall eqs. (D.20) and (D.22) so that, 
on referring to eq. (D.46), 

Ie(Pt^*) = J*(ae"e ''' v tg)\//*dv = ngaglg(i|/*) + \ J V i//*dv 


Similarly, 


Next, 


Similarly, 


= njjaglg(^*) + tpX 


l^(pX) = a^nglgCX) + t^d 


g(ppV/ *) = J (^ q^q + vte)p>/^ *dv 


E. 


= agneIe(pi//*) + te /pdX 




^e"e4 (P’/'*) 


= agnglgCpi//*) + tgpX 


V 

V- 




ixdv 


Ig(ppX) = agfiglgCpX) + tgp0 


-Vgle(X) . 


-teVe(^)- 


(D.63) 


(D.64) 


(D.65) 


(D.66) 


Recalling that eq. (D. 1 2) contains a term multiplied by fl and noting that the term contains 
lg(0), eq. (D. 1 2) involves eq. (D.62); hence, consider 

iJ 2 


= - Ig(R2v/*) ~ Ie(R2) + Ig(X) . 


(D.67) 


Each term of eq. (D.67) is treated separately as follows: 

Ig(R^^*)=^(v^ + ag^ + h^)i//*dv = JvdX + (ag^ + h^) Ig (i//*) 

= vx| 


Jxdv+ (ag2 + h2)lg(i//*) =vX -Ig(X)+ (ag2 + h2)lg(t/^*) 


thus, 




(v3 / 3 ) 


S22lg(0) = 2Ie(X) - (ag2 + h^) - vX 

Similarly, consider 

D.'^Jid) = -h^K^*) + J(X) - ^ aglg(,|/*) 


_+(ag2 + h2)(v+-v-). (D.68) 


(D.69) 


The above results complete the reductions which are possible without becoming heavily 
involved in the theory of special functions. For the treatment of the remaining integrals, 
recourse is made to techniques that involve approximating the integrand with expressions 
that are integrable in closed form. 

D.4 APPROXIMATIONS FOR EVALUATING THE FUNDAMENTAL INTEGRALS 

This section treats in detail the computation of the fundamental integrals le(^), Ie(^)> hJ(\(/), 
J(X), and J(0). The approach which is used is this: Treat singular parts of an integrand 
separately and use a high order polynomial approximation on the regular part of the inte- 
grand (also, functions of R are approximated by polynomials in R). 

For the evaluation of J integrals the operations are as follows: 

The unsteady kernel function, viz., 

j g-iJ2R 


is expressed as 




where 


is the steady kernel function; hence, 


1 

47tR 


= \Ij* = - 1 ) / R. 


(D.70) 


(D.71) 


The function * 1 // is an analytic function of R; consequently, it can be estimated quite 
accurately by polynomial interpolation. 
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In general, the approximation procedure consists of determining the coefficients for the 
approximation having the following form: 

N 

E Cj^ (*'/') R"- 

n = 0 

All that is required of this approximation is that it be accurate for values of R ranging 
over the boundary of the panel. (The reason for this is shown below.) Recalling the 
definition of J, i.e., eq. (D.45), 


*i|/d0-27rCo*)//(lhl); 


since *\}/ accurately approximates *i// on the boundary, 

^ y *V/d0- 27 tCq*i// (I hl) 

{*i//(lhl)-*^^(lhl)} 

N 

n=0 ' ' 

where the presence of the second term allows I) to be a relatively poor approxi- 

mation to *^( h ). 

In a similar fashion, X and d are approximated by the formulae: 

N 

X « X - E Cj^(X)R*^ 

n=0 

N 

d c„( 0 )R" 

n=0 

where a higher degree polynomial is used for 6 than is used for either *\jj or X for 
two reasons: (i) it involves no greater computational cost and (ii) the limiting 
behavior of 0 as 12 tends to zero is given by 

lim0 = R3/3. 


= E f 

® /" 

9t 


(D.72) 


(D.73) 


(D.74) 


(D.75) 


(D.76) 
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thus, if N - 1 , N + 2 = 3 will be required to follow the behavior of d in the limit of small 
frequency. Using the approximations shown by eqs. (D. 75) and (D.76), 

N 

J(X)«^ £ Cj^(X)J(R") -27 tCo {X(lhl)-Xdhl)} (D.77) 

n=0 

N 

J(0)« X) Cj^(0)J(R")-27rCo{0(lhl)-0(lhl)} (D.78) 

n=0 


The task of evaluating the integrals has now been reduced to one of computing hJ(^Q) 
and J(R*^) for n = 0,,..,N+2. Those operations are facilitated by a recursion formula for 
J(R”). Using the definition of J(f), i.e., eq. (D.45), and some of the identities derived in 
the preceding (viz., eqs. (D.l), (D.27), and (D.46)), 

J(r") = X) y* R"dv?-27rColhl" = J (p2 + h2)R«-2 d<p - 2irh2 C^lhP^ 

e 0e" e 0^' 

. + v"^ 

= ^ f p2R"-2- ^ + h2j(Rn-2)^ j agR"-2dv+--k2j(Rn-2) 

e 0e" P e v“ 

r (D.79) 

= E aele(R"“^) + h2j(R"-2). 
e 

For n = 1, note that 

’ (D.80) 


consequently, eqs. (D.79) and (D.80) provide a recursive relation for evaluating J(R*^) for all 
n if J(i//q) and J(l) are known. From eq. (D.46) 


J(l) = 



e 0e“ 


-0 


(D.81) 


because the panel edges form a close curve; thus, only the integral hlft/z^) is required. 

The integral J(i//q) has been studied extensively and these studies (ref 17, appendix J.7, J.8; ref 19 
appendix G) have led to two basic forms, called the standard rationalization and the special rationali- 
zation. Tlie standard rationalization is given (the quantity J of equation J.8. 98, ref. 17 is to be iden- 
tified with hJ(i//Q) j, 

hJ(i//Q) = -sgn(h) j^27T-2(7r-Ph (X(., Y^.)^J (D. 82) 

panel 
corners, c 
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where 


* ^C-1 

Yc= lhlRctc_lXtc 


(D.83) 


and the phase function, Ph(.,.), defined by 

Ph(X^,, Y^) = arg (X^ + iYc),-7r <e<ir, 

is essentially equivalent to the standard intrinsic Fortran function ATAN2. (It is pertinent to remark 
that the expression [ -J(^q) ] of this document is equal to the expression H(l,l,3) that figures so 
prominently in appendixes D and G of ref. 19.) 


The edge integrals, Ig, are determined using eq. (D.79) for values of n=l ,2,...N+2. 

The procedures for Ig(i//) and Ig(X) follow by combining eqs. (D.71) and (D. 72) to obtain 

N 

Cj^(*i//)R". (D.84) 

n=0 

Applying the operator Ig to this relation and to eq. (D.75) yields the following estimates: 

N 

le(>//)^Ie(l/R)+ E Cn(*\I/)Ie(R") (D.85) 

n=0 


leTO- E C„(X)I^(R") ; (D 85, 

n=0 

consequently, the computation required consists of the evaluation of Ig(R^) for 
K = -1,0,.. .,N. This is exactly the same set of information which is required for the 
evaluation of the J integrals, cf. eq. (D.79). 


The following provides a recursion formula for the required integrals: 

v^ Vg 

= R^ds'= /'®p + g^2jRk-2dv 

''e V- • 

for 
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= f dv + gg\ (rJ^-2) 


"e 

/ 


/Rk\ve^ 7 ? / \ 

(it) R^dv + ge2(Rk-2) 


'"e .Ve 


R^ 




Rearranging terms yields the recursion formula as 


where 


I 

( 

v+ 

II 


" +ge^Ie(R‘‘-2) ■ 

\ 

1 


1 

k + 1 

a(vR^^) 



a(vR^)= vR*^ 


,+ 


This recursion relation reduces the requirement to that of evaluating Ig(l/R) and Ig(l). 
These integrals are evaluated as 


Ul) 


7 o‘=' = 


AV 


where 


while 


AV= V+-V; 


Ul/R)« 


r dv _ 

r dv 

R ~ J 

e 

f 

, V 

^ 2.2 

v^+g| 


2 V V -R 


(D.87) 


(D.88) 


(D.89) 
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D.5 POLYNOMIAL ESTIMATION PROCEDURES 

The preceding provides the procedures for evaluating the near field panel influence coeffi- 
cients with the exception of the polynomial estimation procedures for Those 

polynomial estimation procedures are derived in this section. One of the principal consi- 
derations in constructing the polynomial estimates is that they be analytically precise 
in the limit as tends to zero. Since the limiting forms of these functions are 

*V^«(-iJ2)+^(ifi)2R^0 

1 9 

2 

6 ^ R^/3 ^ R^/3 , 

motivation is provided for requiring (at least in the ease of small frequencies) 

-f- = 0(R) 

x-r = o(r^) 

0 -r3/3 =o(r4). 

In fact, it is convenient to require a somewhat stronger condition for *i//, viz., 


+ in (iO)2R = O (r 2) as R 


0 . 


Consequently, if the estimates of *i//, X, and 0 are required for values of R such that 
nR< 1, the estimates can be chosen to have the following forms: 

*5^ - -in n2R + r2pj^_2(R) 

X- R + r2Pj^_2(R) 

0= r3/3 + r4Pn-2(R) 

where Rj^in minimum value of R for which the approximafions mu^ be accurate 

and Pn_ 2 denotes some polynomial (different for each function *i//,X, or0) of degree 


for i2R^in<l 


(D.90) 
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less than or equal to (N-2). Since a polynomial of degree (N-2) is determined by its value 
at (N-1) points, the functions and 6 must be sampled at (N-1) points in order to 

determine the precise forms of the approximations. 


The reason now becomes clear why it is convenient to fit the first two terms of the Taylor 
expansion of and the first four terms of the expansion of X while requiring the degree 
of d to be (N+2) rather than N. With this choice there is a fit to the analytical behavior 
of these kernels for small R while the number of the remaining degrees of freedom is 
identical for each of the three functions. Since the functions of *)//,X, and 9 are all 
structurally similar, all three functions can be evaluated at the same (N-1) points. This mode 
of evaluation requires substantially less computational effort than would be required were 
each function evaluated at a different set of points. 


Eq.(D.90) describes the approximation procedure for the case when ^^R^iin^ 
More generally, the following procedure is used. 


(1) Find the minimum value, R^^iin) and the maximum value, which R 


assumes on the boundary of the panel and define 


max’ 


Ro = 


^min ^^min ^ ^ 


lO ifnR^i„<e 


(D.91) 


where, nominally, e = 0.1. 

(2) The functions *\jj,X,9 are fitted, for R e IR^, , by the polynomial 

forms 


*^ ^*'1' (%) + *'1^ ' K) (R - %)+ (R - Rq) ^ Pn-2 

X = X (R^) + X' (R J (R - R^) + (R . 2 (X) 

j=0 


(3) The remaining (N-1) degrees of freedom are determined by sampling 

X, and 9 at “constrained Tchebysheff absicissae” points chosen to globally 

minimize the expected error in the estimation of and d. The sample 
points, denoted as Rj have the form 
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(D.95) 


' Ro + (Rmax - Ro) f" i = 1 .-.N -I 

where are normalized interpolation points which must be determined by 
solving J nonlinear equations in the following. When the functions have been 
sampled at these points, the remaining degrees of freedom are constrained by 
interpolation. 

The value of N, the degree of the approximation, remains to be chosen by a careful 
error analysis of the above procedure, see sec. D.6. 


Consider the following class of polynomials which are called “constrained Tchebycheff 
polynomials” and which are associated with the optimal choice of interpolation points 
for the constrained interpolation problem, 

f(x) = f(0) + xf'(O) + X'^Pjq_ 2 (x) for 0<x<l. 

A set of basis polynomials is characterized as follows 

Pq = 1, Pj = x,P2 = x2; 

and, for m> 0, Pj ^^.^2 is characterized by the conditions: 


Pm+2<0) = ‘’m+2(0) = Pm+2<» = 


Pm+2 m extrema between 0 and 1 at Xj fori=l,...,m and these extrema satisfy 


Pm+2 


(Xi) = 0 


Pm+2(^i) = (-l)^-^ 

Diagramatically, the basis polynomials appear as 


(D.96) 

(D.97) 

(D.98-a) 

(D.98-b) 

(D.98-C) 
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The problem is to find the extremum points of Pj ^+2 for m and, in this way, charac- 
terize the family of basis polynomials. Once the basis polynomials are known, their roots 
can be determined to obtain the natural analogues of Tchebycheff abscissae. The roots of 
f*n+l ('''hich are N— 1 in number in the interval 0<x<l) provide the evaluation points of 
f(x) for which the optimal solution to the problem is obtained. 


The nonlinear algebraic equation which must be satisfied by the set of extermum points 
for Pjn+2 derived as follows. Using the conditions eqs. (D.98a) and (D.98c), it is clear 
the Piy^+ 2(^) iriay be expressed in terms of the numbers xj by the following formula (Lagrange 
interpolation): 

m 

Pm+2W = E Sj(x)/8j(Xi) + (-l)m5j^^j(x)/ 6j^+j(l) (D.99) 

i=l 


where 


8|(x) = x2(x-l) n(x-xj) 


'm+1 


j- 1 


(x) = x^ n ^x - Xj) . 


j = 1 


(D.lOO-a) 
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The actual positions of the points xj are then determined by improving the constraint shown 

as eq. (D. 98b). Computing Pj^^ 2 (^i)> 


m 


1 k 1 


i=l 

i=^k 


+ (-l) 


k-1 


m 

2 1 .e- 

— + + V 

%-l ^ 


1 


^k 


jzzj % ■ ^i 

iv^k 




(D.lOl) 


where 


^k ^k (^k) ’^m+1 ^m+1*-^^ 

Eq. (D.lOl) motivates the definition of a function F as follows 


(D.102) 


m 


Fk(x) = E 


i=l (%-^i) 
is^k 


(-1)'-* , (-1) 


k-1 


8; 


+ 


5k 


(-1) . (-1) 


k-n 


L5m+1 5j, 


+ (-l) 


k-1 




m+1 


= (-l) 


,k-l. 


2 1 /(-l)i-l _ (-l)k-l' 


+ E 


%5k xj,-Xj 

iv^k 


5; 


m+1 

P(x)=E (- 1 )'^ 5j(x)/6i(xj) 

i=l 


m+1 

6j(x) = x2 .n^ (x-xj) 

j ^ i 
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Syx) _ 2 

6j(x) X 

S'i (^k) 

P'(Xk) 


^m+1 = 1 


. 5|(x) 5i(x) 

p'(x)=E (-1)^-1 


i=l 


6i(x) 8j(xj) 


m+1 


+ E 


•-1 

J=1 J 

j^i 


X' o 1 

^ i(^i) - — +E 
SiC'i) j=i ^i-^j 

m 




. , S|,/6, 

E(-l)‘->/ ‘, + (-1) 

i=l 

k-1 

2 

+ E 


% 

j=i 

i¥^k 


- 

j#k 


rC-l)*"'* 

. (-1)^ 

%^k Xk-Xi 1 

5 

L. 

i 



Xk-Xj 


i^k 
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D.6 ERROR ANALYSIS FOR CONSTRAINED POLYNOMIAL FITTING 

In this section a rigorous error analysis is presented for the polynomial approximations to 
the functions *t//(R), X(R) and 0 (R). The error bounds are relative error bounds; 
for *\p, the error bound takes the form 

H£re, *i//(R) is the transcendental function to be approximated by a polynomial in R, 
*i//(R) is the associated polynomial approximation, and is a “gauge function” 

chosen so that the following approximate bound holds: 

The small positive number 17 measures the relative accuracy of the approximation. It turns 
out that because the function *\jj(R) is analytic in R with all of its derivations nicely 
bounded, a value of r] - 10"^^ can be achieved for all practical purposes by polynomials 
of degree less than or equal to 8. 


Now, in the evaluation of aerodynamic influence coefficients, integrals of the form 


J'J'*\jj ds' and J '■ 


\pd^ 


are approximated by computing the corresponding integrals for the polynomial approxima- 
tions, that is. 


JJ'*\pds' and J' *\pd^ 


Using the relative error bound for *\p given above, the quality of these approximations is 
assessed by means of the following integral error bounds; 


ff *^ds' - jj*i]/ds' 

< 


tjAs' and! j*^dip 

< 

*'P,o 




1 Ej Ee 
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From this, it is easily seen that tj provides an appropriate estimate of the relative accuracy 
of^the integrals as well. Thus, a relative error bound for the polynomial approximation 
*\jj translates directly into a similar relative error bound for the required integrals. 


Having described (in a general sort of way) the nature and the purpose of the error bounds, 
it is now appropriate to present and derive the precise error bounds that are required. As 
a beginning, this presentation and derivation requires the development of some notation. 


The functions *\jj(R), X(R) and 0(R) for which polynomial approximations are desired 
are given by 



(D.103) 


R 

X(R) = J i|/(R')R'dR' 
0 



(-i^2) 


(D.104) 


R 

0(R) = Jx(R')R'dR' = r3 j - 1 - (-iJ2R)] / (-iJ2R)2 

0 

- [e’*^^ - 1 - (-ifiR) -i(-i^2R)2] / (-iflR)3 j . 


(D.105) 


Note that for small frequencies (J2^0), the functions X and 9 exhibit the following 
limiting behavior: 

= -il2 + 0 (D.106) 

X=R + 0(fi) (D.107) 

0 = R3/3 + O(fi) (D.108) 

These limiting properties motivate the definition of the “gauge functions”: X^, 

which will be used later in the statement of the error bounds. The definitions of the gauge 
functions are as follows: 



(D.109) 

X 

o 

ill 

(D.llO) 

0o = R^/3. 

(D.lll) 
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Now, when approximations are desired for the *\p, X, and 0 for R lying in the internal 
[I^o> l^maxl > the approximating polynomials are assumed to have the following structure: 

•fCR) = •((- (R„) + (Ro) (R - Ro) + (R + Ro) ^ Rn-2 I'^HR) 

X(R)-X(R„) +X'(R„) (R-R„) + (R-Ro)^Pn-2 KKR) 

3 

0(R)= E 0^’^(Ro)(R-RoV7j! + (R-Ro)'^PN-2[^](^^)• 

j=o 

Here, P]sj_ 2 [f](R) denotes a polynomial of degree N-2 in the variable R, determined by the 
function f. (Note: N is restricted to be greater than or equal to 1 . For N = 1 , N-2 = -1 
and the convention is adopted that P_] [f] (R) = 0.) The coefficients of the polynomials 
^N-2[*’^ 1 ^N-2t^ J determined by requiring that the polynomials 

*i?, X , e agree with the functions *i//, X, 0 at the N-1 values of R, viz., 

Ri = Ro+(Rmax-Ro)Zi'^^‘’ i=l ,2,....N-1 . 

The numbers i = 1,...,N-1 are the non-zero roots of a polynomial Pjq+i^^) having 

the form 

N-1 

Pn+i(X)=Kn+i n (x-z/n+d), 

i=l 

The polynomial Pj^+|(X) is a “constrained Tschebycheff basis polynomial” which is a 
family of polynomials that will be described below. 


The accuracy of the polynomial approximations, i.e., eqs. (D. 1 12), (D.l 13), and (D.114), is 
determined by two parameters : 

(i) N, which determines the degree of the various polynomial approximations, 
and 

(ii) A0, the total phase variation experienced by the function for 

Rc[Ro)Rmax^ 


where it should be noted the A(^ is given explicitly by 


^^=|^KRmax-Ro)- 


(D.l 12) 
(D.l 13) 

(D.114) 

(D.l 15) 
(D.l 16) 


(D.l 17) 
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For R lying in the interval [Ro’^^max^’ the accuracy of the polynomial approxima- 
tions, i.e., Eqs. (D.l 12), (D. 1 13), and (D.l 14), can be described precisely by the following 
relative error bounds: 

c 

|(*1 ^ -*f)/ *^,o| < (D.118) 
|(X-X)/Xo |<Dn+i(A^)N (D.119) 


1 / , 

3 

|(0-0)/0o 

N + 3 _ 


'"N+1 .Nf+1 

■^%+i 


(A0) 


,N 


The parameters j , Djq+j (tabulated in Table 8) are called “truncation error 
coefficients”; they decrease rapidly with increasing N. The bounds, viz., eqs. (D.l 18) 
(D.l 19), and (D.l 20), clearly show that if N is chosen large enough that 


(D.l 20) 


%+! 

N + 2 


(A0)N+1 


<r? 


(D.121) 


and 

Dn +1 (A0)N<r?, (D.l 22) 

then the following relative error bounds will immediately follow: 

l(X-X)/Xo K?? 
l(0-0)/0o k|t? 


Thus, for any specified error tolerance, 17, an appropriate value of N is selected so 
that the two inequalities: eqs. (D.l 21) and (D.l 22), are enforced, consequently, assuring 
sufficient accuracy is obtained in the polynomial approximations: eqs. (D.l 12), (D.l 13), (D.l 14). 
It is pertinent to remark that for A0 = 2ir / 6 (so that [Ro^Rmaxl spans 1/6 of a wave 
length) and t? = 1 0'^ ^ ( 1 0 significant digits), the inequalities eqs. (D. 1 21 ) and (D. 1 22) are 
satisfied provided N > 8. A fairly complete table of the maximum phase variation, 

A0, allowable for a given N and 77 is presented in Table 9 . 
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Table 8. — Truncation Coefficients * and Zeros ** of Constrained Tschebycheff Polynomials * * * 


N 

2 

3 

4 

5 

6 

C(N+1) 

D(N+1) 

.176488E-01 

.333095E-01 

.78891 5E-03 
.261056E-02 

.327165E-04 

.168479E-03 

.120822E005 

.894424E-05 

.396476E-07 

.399135E-06 

1. 

2. 

3. 

4. 

5. 

.894107456974979 

.606717423814530 

.951856618566651 

.418254069443174 

.761423683974629 

.971712927869628 

.301674740674210 

.590322089237812 

.838296069702794 

.981195629247498 

.226595419540605 

.462063955210008 

.694739789581133 

.882789721841533 

.986532130357020 

N 

7 

8 

9 

10 

11 

O O 
?? 

.116355E-08 

.152913E-07 

.307903E-10 

.511954E-09 

.740598E-12 

.151990E-10 

.163101E-13 
.404951 E-1 2 

.331018E-15 

.997958E-14 

1. 

2. 

3. 

4. 

5. 

£ 6. 

7. 

8. 
9. 

10. 

.175946932336543 

.368204640983135 

.574775757973370 

.765287083020215 

.911009103656216 

.989854498274183 

.140351289637167 

.298873937094264 

.477851060077372 

.656524898416237 

.813554821066919 

.930077552626063 

.992071284053473 

.114455639954652 

.246743790067857 

.401168084345066 

.563272199092783 

.717353855456647 

.848391222077275 

.943584280361694 

.993627507931414 

.095862800941203 

.095062800941203 

.340339265992784 

.485427407681714 

.630434630202327 

.763662228532745 

.874338729060653 

.953508620868845 

.994763511480262 

.080180344720937 

.080180344720937 

.291682130160915 

.420942838028783 

.554690530009402 

.683850675866850 

.799639161635557 

.894174252256846 

.961018408761163 

.995618902517403 


*C(N+1) and D(N+1) 


**L (N+l)!'^''' 

pK 1 2=1 

***Pm+i(X),N=2 11 
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Table 9, — Maximum Admissible Phase Variation * 








N 

4 

6 

8 

10 

12 

2 

.0752 

.0075 

.0008 

.00008 

.0000075 

3 

.4992 

.1075 

.0232 

.0050 

.0011 

4 

1.3244 

.4188 

.1324 

.0419 

,0132 

5 

2.4179 

.9626 

.3832 

.1526 

.0607 

6 

3.6825 

1.7093 

.7934 

.3682 

.1709 

7 

5.0646 

2.6232 

1.3587 

.7037 

.3645 

8 

6.5074 • 

3.6594 

2.0578 

1.1572 

.6507 

9 

7.9904 

4.7902 

2.8717 

1.7215 

1.0320 

10 

9.5116 

6.0014 

3.7866 

2.3892 

1.5075 

11 

11.0543 

7.2730 

4.7852 

3.1483 

2.0714 


^Maximum phase variation, for a polynomial approximation of a given order N 

and relative accuracy 17. 


- log -|Q (97) [77] (significant digits) 



Having stated and discussed the fundamental error bounds : eqs. (D. 118), (D. 1 1 9), (D, 1 20), it 
is now appropriate to present a derivation and proof of their validity. 

As a starting point for this discussion, consider a family of polynomials {Pj } j=o> called 
“Constrained Tschebycheff basis polynomials.” These polynomials are characterized by 
two properties: 

(i) For j > 2, Pj(x) = O(x^) as X ->■ 0 and Pj(x) has j-2 zeroes in the interval 
(0,1) (Constraint property). 

(ii) Pj(x) oscillates between -1 and +1 forxe[0,l] (Tschebycheff property). 


> (These are defined in this way by convention.) 


The first of these two properties, the constraint property, allows one to write down 
immediately the formulae for Pj(x), 

•■0 =' ' 

P,=x 

P2 = x2 ) 

j-2 

Pj=Kjx2n (x- z^) . 

C=1 

The Tschebycheff property now determines the placement of the zeroes zg^^ and the 
value for the normalization constant Kj. In practice, the polynomial Pj(x) is found by 
computing the “external ordinates” : f • 1 j “ 1 

\ i n=\, 

viz., points satisfying the conditions 

xf < . . ■ < Xj_2^> < Xj_/j) = 1 


Pj' (xg^)^ =0 fi = 1 ,2,. . . j-2 (extremal property) 

Pj (xg^^^ = (-1)®'^ £ = 1,2,. . .J-2 (oscillation property). 


(D.123) 


(D.124) 

(D.125) 

(D.126) 


2 

Using the characterization shown by eq. (D.126) together with the constraint Pj(x) = Cx 
for small x, Pj(x) can be written explicitly in terms of 



j-1 
e = i 


by using Lagrange interpolation basis functions (c.f. ref. 19, pp 284-285). One finds that 

j-1 

PjW = E (-1)^'^ (D.127-a) 

£=1 
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where 5g(x) are given, viz., 

j=l 

5g(x) = n (x - Xj.^^). 
r=l 

This representation, coupled with the external conditions eq. (D.l 25), then provide (j-2) 
unknown external coordinates xg^\ £=1 ,... J-2. Solving these equations numerically 
by Newton’s me thod,eq.(D.127)provides an explicit representation of the function Pj(x). 
The representation shown by eq. (D. 1 23), in terms of the zeroes zg^\ can then be 
obtained by solving the equations: 


Pj(z8a))=o 

This has been done and the zeroes C = 1,...,N-1, N = 2,...,1 1 are tabulated in 

table 8. Having computed the condition eq. (D.l 26) with C = j-1 provides an 
explicit and useful formula for the normalization constant Kj. One obtains 

j-2 

(-i)i-2 = Pj(xj9;) = Pj(i) = Kj n (i-zjfl)) 

«=1 

SO that 

j-2 

Kj=(-i)jjn (i-zj<J)) 

£=1 


Having described the constrained Tschebycheff polynomials in sufficient detail, the error 
analysis can now begin. 


Suppose that a function f(x) is approximated on [0,1 ] by a polynomial %(x) of degree 
N by requiring that fjq(x) interpolate the following (N+1) pieces of data: 

f(o), f'(o),f (zg(N+l)^ ,C-1,2,...,N-1. 

Thus, fjq(x) satisfies, 

%(o) = f(0) 
fN>(o) = f'(0) 


(D.127-b) 


(D.l 28) 
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When these conditions are satisfied, the standard error analysis for polynomial interpolation 
then provides the estimates (c.f. ref. 19, p. 278) 


f(x) - f^(^) ~ 


f<N+l)(g) 

(N+D! 




Using eq. (D. 1 23) with j = N+I , one can then write 


f(x) - %(x) = 


f<N+l)( ^) P^+i(x) 
(N+1)! Kn+1 ■ 


Now, since |Pj(x)|<l for all xe [0,1 ], the error bound 

|f(x) - %(x) |<Cj, 4 +l lf<N+l) (^)l, ^e[0,l ] 
follows directly provided the truncation coefficient is defined by 

= 1 / [(N+ 1)! %+l] . 

The bound eq. (D. 1 30) is of fundamental importance in the subsequent error analysis. 
It shows that the error associated with the specified polynomial interpolation process 
is determined by a constant , depending on the order of the fit, and the modulus 
of the (N+l)st derivative of f. The coefficients Cjs^+j decreas^ery rapidly with 
increasing N, an approximate formula being given by 17.1 5x2"^‘®° = 17.15xlO'^' 

Thus, each extra term in the polynomial approximation provides an improvement in 
accuracy of about 1-1/2 decimal digits. 


Another useful bound can be obtained from the estimate eq. (D. 1 29). This bound 
reads as follows 


|f(x)-%(x)|< lx I |f(N+l)(j)| 

where Djq+| is defined by 


_ 1 max 

^N+1 “(N+i)! 0<x< 




(D.129) 


(D.130) 


(D.131) 


(D.132) 


(D.133) 


180 



The bound eq. (D. 132) will be useful in obtaining the relative error bounds for X, viz., 
eq. (D. 1 1 9). The coefficients j also decrease rapidly with increasing N, an approxi- 
mate formula being given by 

Nn+1 = 27.9x 2-4-44 = 27.9x10-1 -34 

Finally, before applying our results to the functions *i//, X and d, one final estimate is 
stated. Suppose that a function g(x) is approximated on [0,1] by a polynomial gisj+2(^) 
of degree N+1 by requiring that gj^+2(^) interpolate the following data: 

g(o), g'(o), g"(o), g"'(o), S = 1,2,..., N-1. 

Then the standard error analysis for polynomial interpolation gives 

g(x) - gjq+2 ~ (x-zl^^*^)-(^-^N-l / > ^e[0,l] (D.134) 

(N+3)! ^ 

Having obtained the fundamental estimates, viz., eqs. (D.l 29), (D.132), (D.134), the error 
bounds, i.e., eqs. (D. 1 1 8), (D. 1 1 9), (D. 1 20), can now be obtained. Consider first the bound 
(D.l 18). Defining the variation of R , i.e., A R, as 


a function R(x) is defined for xe [0,1] by writing 

R(x) = (R„ + ARx) e[R„, J. (D.136) 

Applying the estimate shown by eq. (D. 1 29) to the function f(x), i.e., 

fix) = (R(x)) = - l) / R = (-iS2) [(e-i^^R(x) _ i) / (-ifiR)] (D.137) 


one obtains 


*i//(R(x))-*^^(R(x))| = |f(x)-%(x)| <Cjq+i lfl^l4+l)(|)| for?e[0,l] . (D.138) 

Now )(x) can be easily computed from eq. (D.137) to find 

^N+l)(x) = Hn)I-inAR]N+' (D. 139 ) 
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By a fairly straightforward computation, the identity 




is easily proved. Applying this identity with p = N+1 then yields 


/ \N+1 


-ifiR 


z - i^2R (_i(QR)N+2 


f rN+igfdf = J 


rN+l 


for f = -ifiRr. 

The integral on the right has the following bound : 


rN+l 


< 


j>+l |,-iS2Rr| a, < 


Substituting this result into eq. (D.I39) then yields the bound for viz.. 


|f^N+l)(x)| <|J2llT2ARl^ 


+ 1 1 


N+2 


Substituting this bound for the (N+1 derivative of f into eq. (D. 138) then yields 


< Cn+i IT2| I12AR1 N+1 


N+2 


Recalling that ~ “ IJ2ARI, the relative error bound follows directly as: 




^ N+1 . 

< I A0 

N+2 


|N+1 


as asserted by eq. (D.1 18). 


Next, the error bound, shown as eq. (D. 1 1 9) for X, is proved. Applying the bound 
shown as eq. (D. 132) to the function f(x) defined by 

f(x) = X(R(x)) 


(D.140) 


(D.141) 


(D.142) 


(D.143) 


(D.144) 


(D.145) 
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one obtains 


IX-XKDN+l Xlf<^+^)(d . 


(D.146) 


Now the (N+1)®^ derivative of this new f(x) is given as 


t<N+l)(,)=(i.) '^'"‘ [(e-iS2RW . 1 ) / (-ifl)] = [ARe-MR(x)] 


SO that 


= (AR) (-iJ2AR)^ 


|f(N+l)(x)| <aR(A</))^. 

Substituting this bound into eq. (D.146) and dividing by the gauge function, given by 
Xq = R, yields the following bound: 

|(X-X)/X„|< i Dm+,XAR(A«N = Dn,.,(A«N^^'). 

This result directly implies the bound eq. (D. 1 1 9) by virtue of the inequality 


(D.147) 


xAR _ ' ^o 

R(x) “ R(x) 


= 1 - Rq / R(x) < 1 


(D.148) 


which implies 


|xAR/Rl< 


Finally, the error bound, shown as eq. (D. 1 20), is proved. Applying the estimate eq. (D. 1 34) 
to the function g(x), 

g(x) = 0(R(x)) (D.149) 


one obtains 


N-1 


0(R(x))-0 (R(x)) = - 


(N+3)! 


n (i 




(N+3) 


[0(R(x))l 


C=1 


x=|. 


(D.150) 
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A straightforward computation yields the following equation for the derivative expression: 


,N+3 


N+2 


(dx) [R(x)X(R(x))]AR 


/d\^ 

\ r n /■ \\i / '^ \ 

i 

[x(R(x)) + Re-i^R ] (AR)^ 


d\N+l 


.N 


- (— ) [(2-inR)e'^^^] (AR)3 
= (N + 2 - i^2R)e'i^^ (-iJ2AR)^ (AR)^. 

Substituting this into eq. (D.150) and dividing by the gauge function 6^ = R^/3 yields 


4 N-1 




( ^ - ifiR) e~‘^^ (-iA0)N AR- 

L 1 2 


Breaking this up into two terms as indicated yields 


6 -e _ 3e“*^^(-iA0)N 
~ (N + 3)! 


(N + 2) 


/xAR' 




Zo^N+l) 


+ (-i^2AR)^^^ ^x^ V (x-z/N+1))^ 


Bounding this expression, one obtains 


d-d 


< 


3 • 1 • A0 
(N+3)! 


N 


[(N + 2) • 1 • [(N + 1)! Dn+i] + lA0l • 1 • [(N + 1)! Cn+,]] 


In obtaining this bound, the definition of shown by eq. (D.133) has been used 

together with the bound eq. (D.148) and the following identity: 


max 

X 


x2 n (x-zp(N+l)'j 
C=1 ' ' 


= max 
x 


p (x) 
^N+1 


K 


N+1 


K 


N+1 


(N+1)!Cn+i. 


1 


2 


(D.151) 
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Rearranging the bound eq. (D. 1 5 1 ) slightly now yields the required form of the bound, 
cf. eq. (D.120): 


e -d 


< 


N + 3 


^N+1 


N 


— (A«N+1 

(N+2) 


D.7 FAR-FIELD APPROXIMATIONS 

Recall that the panel influences, viz., eqs. (C.128) through (C.139), describe the potential and 
the gradient of the potential induced at the field point P as the result of source and doublet 
distributions on quadrilateral surface panels. If the panel influence at the field point is 
evaluated, letting the field point recede to greater and greater distances from the panel, 
then the panel influence is found to depend less and less on the details of the panel source 
and doublet distributions; in fact, if the field point is sufficiently far from the panel, the 
panel influence appears to be that of a point source and a point doublet. This observation 
is the fundamental idea underlying the far field approximations to the panel influences. 

The kernel function, \p*, and its gradient, v '/'*) are approximated by a truncated Taylor 
series expansion, and the panel influences are manipulated and separated into the sum of 
a product of two kinds of components. These two types of components are called kernel 
moments and panel moments. The kernel moments depend only upon the field point location 
relative to the point of expansion and are, in particular, independent of the geometry of 
the panel. The panel moments, on the other hand, depend only upon the geometry of the 
panel and may be computed once and for all before any panel influences are computed. 

Because the kernel function , i/' * , exhibits singular behavior as the field point of evaluation 
approaches the point of source or doublet definition, it is clear that the concept for far 
field expansions must break down as the field point approaches the point of kernel function 
definition. As a consequence, the far field expansion is applied only when the distance R 
from the field point P to every point Q on the panel is bounded substantially away from 
zero. Let q^ denote the position of the point Qq lying more or less in the panel center 
and let be the expansion point. Letting p denote the position of the field point, the 
condition for bounding the field point away from the panel is expressed as 

R(p,qo) = I (p-Qo) I >kd 

where D denotes the panel diameter measured in the scaled coordinate system, c.f. sec. 

C. 1. The constant K has been chosen as follows: 

( 1 .5 quadrupole far-fields 
4 dipole far-fields 
12 monopole far-fields 


(D.152) 


(D.153) 
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where the terms monopole, dipole, and quadripole refer to moments of the kernel func- 
tion, and are defined in the following. 


In addition to the condition shown by eq. (C.152), which is adequate to ensure that the function 
1/R is adequately approximated by a truncated Taylor series, it is also necessary to ensure 
that the total phase variation over the panel, viz., 

= (Rj^ax-Rmin)> (D.154) 

is sufficiently small that adequate accuracy is achieved. Roughly, to indicate the nature of 
the requirements imposed by this condition, consider table 10 showing the value of A0 
admissible for the monopole, dipole, and quadrupole approximations when the error in 
the approximation to the function. 




is less than or equal to e. The approximation is given by the truncated Taylor series: 



(D.155) 


and table 1 0 shows the value of A0 for which the approximation attains the accuracy e 
in the interval |^0 q-A 0 , + A0j . Since the absolute value of the error entailed in the 

approximation is provided by Lagrange remainder theorem, whence 



N 


E 

j=o 






1 



(N+l)! 


N+l 


the admissible range for 0 for a given N and a given e is 

A0 = [(N+ 1)! e] _ 


Table 10 is based on eq. (D.l 56). 


(D.156) 


These remarks are aimed at describing the accuracy of the procedure; the following 
describes the details of the procedure itself. 
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Table 10. — Admissible Panel Phase Variation * Related to Truncation Errors 



.001 

.003 

.01 

.03 

.1 

.3 

1, 

monopole 

.045 

.077 

.141 

.245 

.447 

.775 

2, 

dipole 

.182 

.262 

.391 

.565 

.843 

1.216 

3, 

quadrapole 

.394 

.518 

.700 

.921 

1.245 

1.638 


* Value of A0 such that, for (0^ - A(^) < 0 < (0^ + A4>), the approximation 



is accurate to within e 
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As stated above, the underlying principle of the far-field approximation follows from having 
replaced the actual kernel functions by truncated Taylor series approximations; therefore, as 
a starting point for describing the far-field approximation procedure, consider the Taylor 
series expansions for the kernel function and its gradient, viz., 

\l/* = J— 

47tR 


and 


V V*, 


where these expressions are in terms of the scaled coordinates introduced by sec. C. 1 . 
Proceeding in such a way as to preserve maximum symmetry, let 


\p* = \p* (lp-q|)= i//* ^1 p - q'o - (q -%)l) 


where p is the position of the field point P and q is the position of the surface point 
Q. The surface point Qq with position will be the point of expansion; thus, letting 


P = P-qo andQ^q-q^, 


• vl/2 


P = (pigyP-i and Q = (Q‘gijQ> ) ' ' " , 
Pj = g}jPj and Qi = gjjQ> , 

R= [(P'-Qi)gij(pj -0^')]^^^ 


leads to the following expansion of the kernel function: 


Similarly, 


,//*(R) = <//*(P)-Q'^ >//*(P)+-Q‘Q’ — 

api 2 apigpj 


9)//* d\jj* ^ Pk 9i// 


: ^*(P) . 


V i// * = ~ = -a*' ~T 

9Q^ 9P^^ 


9P^ 


(D.157) 


(D.158) 

(D.159) 

(D.160) 

(D.161) 


(D.162) 


= apg- 


£k 


^H=(p)+Qi 

dpi^ 


gpkgpi 2 ap^apiapJ 


i//*(P) 


(D.163) 
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Introducing the following identities: 


9i//* 

3pi 


9P 9i//*_ 9 / g P*^ 9i//* _ 9</^ * 

W ^ ~^™ir 8P ~"P3P 


9pi9pj 

d^yp* 

9pi9pj0P^^ 


j_ 

®ij p 9P 


-~Wk 




2 

xp* 


and defining the operator D by 

and substituting it into eqs. (D. 1 62) and 
kernel function and its gradient: 

xp*m = p*(P) + Q%Dp* +^Q^(^(p^?p^P*-g^p\p*y (D.164) 

v'i//*='agg®^ |pj^Di//*+Qi |^PjPj^D^i//*-gij^D'l'*j 

+ ~QiQ) [PiPjPkD^.^=^ v(gijPk+gikPj+gjkPi)D2>^*]^ (D.165) 

where the following identity has been used : 

Qjpi = QiPi. 

Eqs. (D.164) and (D. 165) lead to the definition of quantities which are termed kernel function 


moments, viz., 

G(o) 

= \p*(P) 



Gb) 

= P^D\p* 



^h) 

= P*pj(D^i//*) -g^Di//* 




= pipjp'^ fD^V'*) - (pV*^ + pjg'*^P'^g‘-’) D^\p* 

(D.166) 


P 9P 

(D. 163) leads to the following truncated expansions of the 


where each of these moments has the property of a tensor which is completely symmetric 
in its indices. In terms of the kernel moments the truncated expansions become 


^W = G(0) + QiG{i) +-QiQjCy 


( 2 ) 


(D.167) 


= % (*^0) + Qi<^(2) + j QiQjGj'l) 


(D.168) 


Having obtained these approximations for the kernel function and its gradient, the panel 
influences shown by eqs. (C.128) through (C. 139) can be evaluated for the far-field approxi- 
mation. The result of that evaluation is expressed in terms of the following panel moments: 


] 

/V 

1 ^ 

Ein 


Qi ^ 

F-- 
ijn / 

Se 

'Wknl 

c c 

^ 1 ' 


- JJ -n 

Qi 

1 

.Wijkn, 

2e 



^n ”n 


nj^ds for ju* =53 mj^iUj^ 


(D.168) 


(D.170) 


'F6nl 


1 1 

Pfiin 


Fsijn 

2e 



ds for o; = 1 ,2 
0 = 1,2 
5 = 1,2 


(D.171) 


where it should be noted that the entries of the F-moments are vectors; morever, all 
these moments are symmetric in the ij indices. 


Using the above definitions, the panel influences are written in the following concise form: 

■'■s' Z /Jsn (G(0)+QiG(|)+5QiQjG("2))dsJ„ 

= E(G(0)En + Gu)Ein + G;i)Eij„)o„ 
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” Se 

x-' / k ik ^ijk \ _ 

= E(G(,)Wt„ + G(2)Wft„ + G(4)Wij|,„j(i„ (D.173) 


C C / k ki 1 kii ^ . 

n 2g 

^G(-|)Ej^ + G(2) Ejj^ + G^4yEj|^j^^ajfj (D.174) 

n 


„/ k ik ijk\ 

VDg =Z^r5„G(j) + Fgi„G(-2) + F§jjj,G(-4)j nj^dsMn 
n ' ' 

''E>3 Z) -(pSn^O ) ■*■ f’5inG(2) f^Sijn ^( 4 )) (D.175) 

n ' ' 
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APPENDIX E 

SECOND ORDER DERIVATIVES OF THE VELOCITY POTENTIAL 


The following terms appears in the boundary condition shown as equation eq. (31): 

(D*-vWs)-fts 

and this term contains the operations: 

V Wg = V - cMq^c - v0s) (E.l) 

which, in turn, contains second order derivatives of the velocity poten tial of the steady 
mean component of flow. The objective of this appendix is to present a means for 
evaluating those derivatives. 


Recall that the doublet strength at a network surface point is equal to the jump in potential 
across the network surface, c.f., eq. (C.27). In the case of a thick body on which Morin o- 
type boundary conditions are applied, c.f., sec. 3.4, the potential vanishes at the interior of 
the body ; hence, the doublet strength at the body surface is equal to the value of the potential at 
the exterior side of the surface of the body, i.e., 

at the body surface. The steady flow panel method computes the doublet strength as a 
quadratic function of the surface coordinates, viz., 

1 .9 • 9 

Ms(Q) = Mo + + 2 2 

hence, the following second order derivatives of the velocity potential are readily evaluated: 


(^s) 

(E.4) 

Of the remaining second order derivatives can be evaluated using the flow equation, i.e., 

V 2<^s-Mo2S- v(c-v0s) = 0. (E.5) 
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Letting 




represent the unit base vectors in the directions of the surface coordinates, it follows that 





The remaining second-order derivatives can be evaluated using the steady flow boundary 
condition applied at aerodynamic surfaces, viz., 

W^.^s = 0, 

which, when expressed in terms of the surface coordinates, becomes 

’ f T T ^ ^ 

0f--UoC-n^. 

The surface is to be represented as a quadratic function of the surface coordinates at each 
panel, i.e., 


where 


and the vector normal to the surface is 

A __ vf 

where 


I vf I =V ^ + (?t7 + ^rjrjV + + 1 • 


(E.6) 


(E.7) 


(E.8) 


(E.9) 


(E.IO) 


(E.ll) 
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Finally, the two remaining second-order derivatives are formed by differentiating eq. E.7 
to obtain the following: 

% = -UoC- [(-f^^%-feA)(w|vf|)-(vf)(|vf|^)/^^^ 


% = - Uo c • [ ( - - ^rivh) (l / 1 ) -(^f ) ( It?)/ ( 1^ f I) 


Eqs. (E.4), (E,6), (E,l 2), and (E.1 3) provide a means for evaluating all nine components 
of the gradient of the mean steady flow velocity potential. 


(E.12) 


(E.13) 
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